Wavelet depulsing of ultrasound echo sequences

ABSTRACT

Methods and associated apparatuses for imaging a target. An echo sequence image of the target is acquired and a log spectrum of at least a portion of the echo sequence image is computed. A point spread function is estimated by one of two methods. According to the first method, a low-resolution wavelet projection of the echo sequence log spectrum is used as an estimate of the log spectrum of the point spread function. According to the second method, an outlier-resistant wavelet transform of the echo sequence log spectrum is effected, followed by soft-thresholding and an inverse wavelet transform. Under both methods, a frequency domain phase of the point spread function also is estimated. the relevant portion of the echo sequence image is deconvolved using the estimated point spread function.

FIELD AND BACKGROUND OF THE INVENTION

[0001] The present invention relates to medical imaging and, moreparticularly, to a method of deconvolving an ultrasonic echo sequence,and to an ultrasound imaging apparatus that employs this method.

[0002] Because of the coherence of the back-scattered echo signals,images obtained from echo ultrasound imaging systems have extremelycomplex patterns that bear no obvious relationship to the macroscopicproperties of the insonified object. The vast majority of biologicaltissues are extremely small on the scale of an acoustic wavelength.Consequently, a signal obtained within a resolution cell consists ofcontributions of many independent scatterers. Interference of thesede-phased echoes gives rise to a pattern that has the appearance of achaotic jumble of “speckles”, known as speckle noise. The specklepattern consists of a multitude of bright spots where the interferenceis highly constructive, dark spots where the interference isdestructive, and brightness levels between these extremes. The presenceof speckle noise in an ultrasound image reduces the ability of a user toresolve fine details. Speckle noise obscures very small structures, forexample, early stage tumors, and decreases the reliability of tissuecharacterization. Therefore, the suppression of speckle noise is animportant component of medical ultrasound imaging.

[0003] For the purpose of modeling the interaction of biological tissuewith ultrasonic waves, the biological tissue is considered to be anassembly of reflectors and scatterers. A reflector is a plane interfacethat is large compared to the wavelength and that reflects portions ofthe transmitted energy back towards the transmitter. A scatterer is anobject that is small compared to the wavelength and that scatters thetransmitted signal in all directions. Such a system often is modeled asa (most generally 3D) function called the spatial response of insonifiedmaterial, the reflectivity function, or (in medical applications) thespatial tissue response.

[0004] An ultrasound radio frequency (RF) image can be considered toconsist of 1D echo sequences, also known as “RF lines”. Assuming thetissue properties to be uniform in the plane perpendicular to thescanning beam, an acquired 2D RF image can be viewed as the result ofthe convolution of the 2D reflectivity function (which accounts forinhomogeneity in the scanning plane) and the 2D transducer point spreadfunction (PSF). Thus, the RF image can be considered to be a distortedversion of the true reflectivity function, where the distorting kernelis the transducer PSF. This distortion includes the speckle noisediscussed above.

[0005] In principle, it should be possible to measure the PSF in acalibration procedure, and then to deconvolve the PSF from the RF image.In practice, however, this is not possible, for several reasons. Perhapsthe most important reason is that the absorption of ultrasound energy intissues increases with frequency. This frequency-dependent attenuationcauses both the PSF amplitude and the PSF shape to depend on depth inthe tissue, leading to the observed non-stationarity of RF sequences.

[0006] In medical ultrasound, a pulse is transmitted into the tissue tobe imaged, and the echoes that are backscattered to the emittingtransducer are detected as a voltage trace RF line. The RF lineconventionally is modeled as being a convolution of a hypothetical 1DPSF with a hypothetical 1D tissue reflectivity function. Assuming thatthe scatterers on each image line are located on a uniform grid and thatthe system impulse response is range shift invariant along each imageline, a discretized version of the received signal can be written as:

rf[n]=a[n]*s[n]+noise[n]  (1)

[0007] where n is a time index, rf[n] is the RF line, s[n] is thetransmitted ultrasound PSF, a[n] is a reflectivity sequencecorresponding to the reflectivity function, noise[n] is measurementnoise, and “*” represents convolution. Because the frequency-dependentattenuation process appears as a decrease with distance of the meanfrequency and amplitude of the PSF, it is commonly assumed that thereceived echo signal may be expressed as a depth-dependent PSF convolvedwith the tissue reflectivity function. To make the PSF “locationdependent”, s[n] in equation (1) is replaced by s[n,k], where k is thelocation index. This leads to the observed non-stationarity of the RFlines received from the tissue. In order to deal with thisnon-stationarity, the RF-sequence is broken up into a number of possiblyoverlapping segments, such that within each segment thefrequency-dependent attenuation process can be ignored and equation (1)holds. The problem of tissue characterization is thus reduced to a setof blind deconvolution problems: for each segment of a given RF line,the respective ultrasonic PSF should be estimated and removed.

[0008] To this end, homomorphic signal processing has been applied torf[n]. Ignoring the noise term on the right hand side of equation (1)for now, transforming equation (1) to the frequency domain gives:

RF(w)=A(w)S(w)  (2)

[0009] i.e., in the frequency domain, the frequency spectrum A(w) of thereflectivity sequence is multiplied by the frequency spectrum S(w) ofthe PSF to give the frequency spectrum RF(w) of the echo sequence. Thesespectra can be written as

RF(w)=|RF(w)|e ^(j·arg(RF(w)))  (3)

A(w)=|A(w)|e ^(j·arg(A(w)))  (4)

S(w)=|S(w)|e ^(j·arg(S(w)))  (5)

[0010] Taking the complex logarithm of both sides of equation (2) thengives

log|RF(w)|=log|A(w)|+log|S(w)|  (6)

arg[RF(w)]=arg[A(w)]+arg[S(w)]  (7)

[0011] As described in the Annexes, the log spectrum of the echosequence, log|RF(w)|, thus is a sum of a smooth and regular logspectrum, log|S(w)|, of the PSF and a jagged and irregular log spectrum,log|A(w)|, of the reflectivity sequence. Cepstrum-based techniques havebeen used to exploit the qualitatively different natures of the PSF andreflectivity log spectra to isolate the PSF log spectrum for the purposeof estimating the PSF and deconvolving the estimated PSF from the echosequence. See, for example, Torfinn Taxt, “Restoration of medicalultrasound images using two-dimensional homomorphic deconvolution”, IEEETransactions on Ultrasonics, Ferroelectrics and Frequency Control, vol.42 no. 4 pp. 543-554 (July 1995); Torfinn Taxt, “Comparison ofcepstrum-based methods for radial blind deconvolution of ultrasoundimages”, IEEE Transactions on Ultrasonics, Ferroelectrics and FrequencyControl, vol. 44 no. 3 pp. 666-674 (May 1997); J. A. Jensen and S.Leeman, “Nonparametric estimation of ultrasound pulses”, IEEETransactions on Biomedical Engineering, vol. 41 pp. 929-936 (1994); andJ. A. Jensen, “Deconvolution of ultrasound images”, Ultrasonic Imaging,vol. 14 pp. 1-15 (1992). Cepstrum-based techniques, however, suffer fromcertain limitations, as discussed in Annex A. Briefly, the complexcepstrum of a signal of finite duration has been shown to extend toinfinity. This invariably leads to aliasing errors when the DiscreteFourier Transform or a similar discrete numerical method is used tocompute the cepstrum.

[0012] There is thus a widely recognized need for, and it would behighly advantageous to have, a method of estimating an ultrasound PSFthat would overcome the disadvantages of presently known methods asdescribed above.

SUMMARY OF THE INVENTION

[0013] The present invention includes two innovative methods forestimating the PSF of an echo sequence. The first method is based on alow-resolution wavelet projection of the log spectrum of the echosequence. The second method is based on the “soft-thresholdingde-noising” algorithm of David L. Donaho and his coworkers, modified toaccount for the fact that the log spectrum of the reflectivity sequenceis not normally distributed.

[0014] According to the present invention there is provided a method ofimaging a target, including the steps of: (a) acquiring an echo sequenceimage of the target; (b) computing a log spectrum of at least a portionof the echo sequence image; (c) computing a low-resolution waveletprojection of the log spectrum; (d) estimating a point spread functionfrom the low-resolution wavelet projection; and (e) deconvolving the atleast portion of the echo sequence image with the point spread function.

[0015] According to the present invention there is provided a method ofimaging a target, including the steps of: (a) acquiring an echo sequenceimage of the target; (b) computing a log spectrum of at least a portionof the echo sequence image; (c) effecting an outlier-resistant wavelettransform of the log spectrum, thereby producing a plurality of waveletcoefficients; (d) soft-thresholding the wavelet coefficients; (e)applying an inverse wavelet transform to the soft-thresholded waveletcoefficients to obtain a point spread function log spectrum; (f)estimating a point spread function from the point spread function logspectrum; and (g) deconvolving the at least portion of the echo sequenceimage with the point spread function.

[0016] According to the present invention there is provided an apparatusfor imaging a target, including: (a) a transducer for acquiring an echosequence image of the target; and (b) a processor for: (i) computing alog spectrum of at least a portion of the echo sequence image, (ii)computing a low-resolution wavelet projection of the log spectrum, (iii)estimating a point spread function from the low-resolution waveletprojection, and (iv) deconvolving the at least portion of the echosequence image with the point spread function.

[0017] According to the present invention there is provided an apparatusfor imaging a target, including: (a) a transducer for acquiring an echosequence image of the target; and (b) a processor for: (i) computing alog spectrum of at least a portion of the echo sequence image, (ii)effecting an outlier-resistant wavelet transform of the at least portionof the log spectrum, thereby producing a plurality of waveletcoefficients, (iii) soft-thresholding the wavelet coefficients, (iv)applying an inverse wavelet transform to the soft-thresholded waveletcoefficients to obtain a point spread function log spectrum, (v)estimating a point spread function from the point spread function logspectrum, and (vi) deconvolving the at least portion of the echosequence image with the point spread function.

[0018] According to both methods of the present invention, an echosequence image of the target is acquired. As understood herein, an echosequence image is a set of RF lines, acquired in parallel, from whichthe final video image of the target is to be calculated. Optionally, theecho sequence image is partitioned into a plurality of segments that maybe either disjoint or overlapping. Subsequent processing is appliedeither to the echo sequence image as a whole or to one or more of thesegments separately. This processing begins with the computation of thelog spectrum of the time series (the whole echo sequence or the segmentthereof) being processed.

[0019] At this point, the two methods of the present invention divergetemporarily. According to the first method, a low-resolution waveletprojection of the log spectrum is computed. Preferably, this waveletprojection is based on a Coiflet wavelet, on a minimum-phase Daubechieswavelet, on a symmetric Daubechies wavelet or on a biorthogonal wavelet.Optionally, this wavelet projection includes an outlier-resistantwavelet transform of the log spectrum.

[0020] According to the second method, an outlier-resistant wavelettransform is essential, not merely optional. Specifically, according tothe second method, an outlier-resistant wavelet transform of the logspectrum is effected, and the resulting wavelet coefficients aresubjected to the “soft-thresholding” algorithm of David L. Donaho etal., or an equivalent algorithm. The outlier-resistant wavelet transformpreferably is based on minimizing an L₁ norm. Most preferably, theoutlier-resistant wavelet transform is the Smoother-Cleaner WaveletTransform of Andrew G. Bruce et al., or an equivalent transform. Robustresiduals may be removed at all resolution levels of theoutlier-resistant wavelet transform or at only some resolution levels.In particular, robust residuals may be removed at only the firstresolution level. An inverse wavelet transform is applied to thesoft-thresholded wavelet coefficients, thereby producing a PSF logspectrum.

[0021] Either the low-resolution wavelet projection of the first methodor the PSF log spectrum of the second method is an estimate oflog|S(w)|. The remaining steps of the two methods are identical. Anestimate of the frequency domain phase arg[S(w)] of the PSF is obtained,preferably under the assumption that the PSF is a minimum phasesequence. Optionally, a Wiener filter is applied to the estimate of|S(w)|. Equation (5) now gives an estimate of S(w), which is inverseFourier transformed to give an estimate of the PSF. Finally, theestimated PSF is deconvolved from the time series (the whole echosequence or the segment thereof) being processed preferably using anapproximate inverse.

[0022] An apparatus of the present invention includes a transducer foracquiring an echo sequence image of the target and a processor forprocessing the acquired echo sequence image according to one of the twomethods.

[0023] Although the present invention is illustrated herein withreference to its primary application to ultrasound imaging, it is to beunderstood that the scope of the present invention extends to anyimaging modality based on receiving echoes of signals transmitted by animpulsive energy source to a target, and to which the principles of thepresent invention are germane (for example, seismic exploration).

BRIEF DESCRIPTION OF THE DRAWINGS

[0024] The invention is herein described, by way of example only, withreference to the accompanying drawings, wherein:

[0025]FIG. 1 is an overall flowchart of the methods of the presentinvention;

[0026]FIG. 2 is a flowchart of the low-resolution wavelet projectionprocedure;

[0027]FIG. 3 is a flowchart of the second method of estimatinglog|S(w)|;

[0028]FIG. 4 is a high level block diagram of an ultrasound imagingapparatus of the present invention.

DESCRIPTION OF THE PREFERRED EMBODIMENTS

[0029] The present invention is of methods, and associated systems, foracquiring an echo sequence image of a target using an impulsive sourceof energy, and then estimating and deconvolving the point spreadfunction of the impulsive source. Specifically, the present inventioncan be used to acquire medical ultrasound images with reduced specklenoise.

[0030] The principles and operation of echo sequence image acquisitionand processing according to the present invention may be betterunderstood with reference to the drawings and the accompanyingdescription.

[0031] Referring now to the drawings, FIG. 1 is an overall flowchart ofthe methods of the present invention.

[0032] In block 10, a 1D ultrasound echo sequence is acquired in theconventional manner.

[0033] In block 12, the acquired echo sequence is optionally partitionedinto two or more subsequences. As described in Annex A with reference toequation (2.1) (and with a slight change in notation), biological tissueacts as a low pass acoustic filter, so that the frequency content of anultrasound pulse decreases with time. To compensate for this effect, theecho sequence rf[n] usually is partitioned into a set of subsequencesr[n,k], such that, in each subsequence, the ultrasound PSF s[n,k] can beconsidered to have a well-defined frequency spectrum. The subsequencesrf[n,k] may be either disjoint (i.e., nonoverlapping), or may overlap toa certain extent. In what follows, the subsequence index k will besuppressed for notational clarity.

[0034] In block 14, the PSF is estimated. The present invention includestwo methods of estimating the PSF. More precisely, the present inventionincludes two methods of estimating log|S(w)|. According to the firstmethod, which is described in detail in Annex A, a low-resolutionwavelet projection of log|RF(w)| is taken to represent log|S(w)|.According to the second method, which is described in detail in Annex B,log|RF(w)| is treated as a noisy representation of log|S(w)|, and thenoise is removed using the “soft-thresholding de-noising” algorithm ofDavid L. Donaho and his coworkers, modified to account for the fact thatthe samples of log|A(w)| are not distributed according to a normal(Gaussian) distribution. Under either method, arg[S(w)] is estimated asdescribed below. Equation (5) then gives an estimate of S(w), whoseinverse Fourier transform is the PSF.

[0035] In block 16, the estimated PSF is deconvolved from the echosequence to give an estimate of the reflectivity sequence, as describedbelow.

[0036]FIG. 2 is a flowchart of the low-resolution wavelet projectionprocedure.

[0037] In block 20, the wavelet used in the wavelet decomposition oflog|RF(w)| is selected. Suitable wavelets include, among others, Coifletwavelets, minimum-phase Daubechies wavelets, symmetric Daubechieswavelets and biorthogonal wavelets, with Coiflet wavelets beingpreferred.

[0038] In block 22, the resolution level of the wavelet transform thatdistinguishes log|A(w)| from log|S(w)| is selected. One way of doingthis is to calibrate the ultrasound probe by recording the pulse emittedby the probe, either directly or after reflection from a single planarreflector, treating this recorded pulse as a proxy for s[n], computingthe corresponding log|S(w)| and determining the number of waveletcoefficients needed to represent this proxy for log|S(w)| with thedesired degree of accuracy. As mentioned above, this measured proxy fors[n] can not be used directly for deconvolution of an acquiredultrasound echo sequence, for several reasons. First, the actual s[n]generally is low-pass filtered by the biological tissue. Second, theactual s[n] always differs from the proxy s[n] because of invariablevariations in the conditions of acquisition, such as the degree ofacoustical impedance mismatch between the ultrasound probe and thetarget. Nevertheless, the proxy of log|S(w)| can be used to define theresolution level of the wavelet transform because the wavelet resolutionlevels needed to accurately represent the actual log|S(w)| and the proxyof log|S(w)| generally are identical.

[0039] The low-resolution wavelet projection itself is defined inequation (3.10) of Annex A. With a slight change of notation, thisequation is: $\begin{matrix}{{\log {{{RF}(w)}}} = {{A_{2^{J}}\log {{{RF}(w)}}} + {\sum\limits_{j \leq J}{D_{2^{J}}\log {{{RF}(w)}}}}}} & (8)\end{matrix}$

[0040] Note that resolution level decreases with increasing resolutionindex j. The first term on the right hand side of equation (8), A₂,log|RF(w)|, is the low-resolution wavelet projection of log|RF(w)| atresolution level J that is computed in block 24. The remainder of theright hand side of equation (8) is the high-resolution portion oflog|RF(w)|. A₂, log|RF(w)| is used in subsequent processing as anestimate of log|S(w)|.

[0041] The second method of estimating log|S(w)| is an extension of the“soft thresholding de-noising” algorithm that was developed by David L.Donaho and his coworkers, and that they have described in a series ofpublications: David L. Donoho, “De-noising by soft-thresholding”,Transactions on Information Theory, vol. 41 no. 3 pp. 613-627 (1995);David L. Donoho and Iain M. Johnstone, “Adapting to unknown smoothnessvia wavelet shrinkage”, Journal of the American Statistical Association,vol. 90 pp. 1200-1224, (1993); David L. Donoho and Iain M. Johnstone,“Ideal de-noising in an orthonormal basis chosen from a library ofbases”, C. R. Acad. Sci. Paris, Series I, vol. 319 pp. 1317-1322 (1994);and David L. Donoho and Ronald R. Coifman, “Translation-invariantde-noising”, Technical Report 475, Department of Statistics, StanfordUniversity, (May 1995). This algorithm, and equivalent algorithms, arereferred to herein collectively as “soft-thresholding” algorithms. Thepurpose of these algorithms is to clean up the noise of a noisy signalin order to extract the signal. In the present context, log|S(w)| istreated as the signal part of log|RF(w)| and log|A(w)| is treated as thenoise part of log|RF(w)|: these algorithms are applied to log|RF(w)| toclean up log|A(w)|, thereby isolating and extracting log|S(w)|. Thesealgorithms assume that the noise obeys a Gaussian distribution. It isshown in Annex B that log|A(w)|, treated as random noise, is expected toobey a non-Gaussian distribution, specifically, a Fisher-Tippetdistribution. To account for this, the wavelet coefficients oflog|RF(w)|, that are input to the soft-thresholding algorithm, arecomputed using an outlier-resistant wavelet transform. Essentially,wavelet projections are used that minimize the L₁ norm instead-of the L₂norm at each level of resolution. In practice, for computationalefficiency, the “Smoother-Cleaner Wavelet Transform” (Andrew G. Bruce,David L. Donoho, Hong-Ye Gao and R. Douglas Martin, SPIE Vol. 2242 pp.325-336 (1994)) is used in preference to brute force L₁-normminimization.

[0042]FIG. 3 is a flowchart of the second method of estimatinglog|S(w)|.

[0043] As in the low-resolution wavelet projection of FIG. 2, the firststep (block 30) is the selection used in the discrete wavelet transformof log|RF(w)|. The outlier-resistant discrete wavelet transform itselfis effected in block 44, as a loop over resolution levels. The loop isinitialized in block 32 at the coarsest resolution level. At eachresolution level, wavelet coefficients are calculated in block 34,robust residuals are calculated in block 36, and the thus-identifiedoutliers are removed in block 38. In block 40, the loop terminationcondition is tested: if the finest desired resolution has not beenreached, the resolution level is made finer (block 42) and the loop isrepeated; otherwise the loop is terminated.

[0044] As discussed above and in Annex B, the robust residualcalculation of block 36 amounts to computing the wavelet coefficients byminimizing an L₁ instead of an L₂ norm at each level of resolution. Thisminimization can be done by brute force. It is preferable, however, fornumerical efficiency, to use the “Smoother-Cleaner Wavelet Transform” ofBruce et al. for this purpose.

[0045] Following the outlier-resistant wavelet transform of block 44,the soft thresholding algorithm of Donaho and his coworkers is appliedto the resulting wavelet coefficients in block 46, and thesoft-thresholded wavelet coefficients are subjected to an inversediscrete wavelet transform in block 48 to give the final estimate oflog|S(w)|.

[0046] Although FIG. 3 shows robust residuals being calculated andoutliers being removed at all the resolution levels considered, it hasbeen found that it is not strictly necessary to remove outliers at allresolution levels. In fact, it often suffices to remove outliers only atthe first (coarsest) resolution level.

[0047] It will be appreciated that the outlier resistant wavelettransform of block 44 also may be used in the low-resolution waveletprojection procedure of FIG. 2, in the projection step of block 24.

[0048] Optionally, following the exponentiation of the estimate oflog|S(w)| obtained by either method, the Wiener filter defined byequation (3.11) of Annex A is used to refine this exponentiated estimateto provide a refined estimate of |S(w)|.

[0049] |S(w)| having been estimated, arg[S(w)] now is estimated. Asdescribed in Annex A, the preferred estimate of arg[S(w)] is a minimumphase estimate. If it is assumed that the PSF is a minimum phasesequence, then log|S(w)| and arg[S(w)| are a Hilbert transform pair, sothat an estimate of arg[S(w)] can be derived from the estimate oflog|S(w)|.

[0050] With both |S(w)| and arg[S(w)] now estimated, equation (5) givesan estimate of S(w). The inverse Fourier transform of this estimate ofS(w) is an estimate of s[n], which is deconvolved from rf[n] byapproximate inverse methods, as described, for example, in A. K. Louis,“Approximative inverse for linear and some nonlinear problems”, InverseProblems, vol. 11 no. 6 pp. 1211-1223 (1995).

[0051]FIG. 4 is a high level block diagram of an ultrasound imagingapparatus 50 of the present invention. Apparatus 50 includes aconventional ultrasound transducer 52 coupled to a digital processor 54via a D/A converter 56 and an A/D converter 58. Processor 54 generatesdigital ultrasound waveforms that are transformed to analog signals byD/A converter 56 and transmitted into the biological target to be imagedby transducer 52. The consequent ultrasound echoes are received bytransducer 52, digitized by A/D converter 58 and processed by processor54 in accordance with the principles of the present invention to provideimages of the target with reduced speckle noise.

[0052] While the invention has been described with respect to a limitednumber of embodiments, it will be appreciated that many variations,modifications and other applications of the invention may be made.

ANNEX A “Nonparametric Estimation of Ultrasonic Pulses Using the FWT andReflectivity Function Acquisition by the Approximate Inverse Technique”Dan Adam and Oleg Michailovich

[0053] 1. Introduction

[0054] In order to describe scattering of the ultrasound from softtissues several theoretical approaches exist. One way is to consider theacoustic inhomogeneity as an assembly of reflectors and scatterers. Theformer is a plane and large, compared to the wavelength, interfacereflecting portions of the transmitted energy back toward thetransmitter, and the latter is a small object, compared to thewavelength, spreading the transmitted signal in all directions,reflecting small portions of the energy. Such a system can be generallymodeled as a some function supplying us for each point of interrogatedtissue reflection strengths, and it is a measure of how much of thetransmitted energy is reflected by an imaginary point reflector locatedat the given point. This function is frequently called the spatialresponse of insonified material and it strongly depends on thetransmitter aperture position and geometry. Further we will refer to itas the spatial tissue response or reflectivity function. With certainassumptions, which will be discussed latter, any receivedradio-frequency (RF) image can be modeled as a convolution product ofthe spatial impulse response of the observed tissue and the transducerpoint spread function (PSF). Note that both the tissue reflectivityfunction and the transducer PSF are generally 3-D functions, but thedimension of the problem can be decreased by several ordinaryassumptions, which might be made without affecting the modelpreciseness.

[0055] In order to improve the transducer resolution, diverse kinds ofapertures (namely, curved, phased arrays) are used. Such devices haveelectronic focusing in the lateral direction (transverse to the beam)and by curving the piezoelectric elements a fixed (not electronic) focusin the elevation plane also can be obtained. The ultrasound beamscanning is, usually, performed in the horizontal (transducer) plane.For such a scheme, the received RF-image G(n,m) can be modeled as a 2-Dconvolution of the 2-D transducer PSF S(n,m) and the 2-D tissuereflectivity function F(n,m)[1-2, 4, 6].

G(n,m)=S(n,m)*F(n,m)+U(n,m)  (1.1)

[0056] Here U(n,m) is the white noise taken from the normal distributionwith zero-mean and variance σ_(u) ² that always has to be added for realsystems; the indices n and m are lateral and depth (axial) increments,respectively. Note that the mapping from 3-D space to 2-D space leads tonon-uniqueness of solution in general. Yet if we assume that all imagefeatures in the imaging plane extend (in limits of an effective heightof the beam) perpendicular to the corresponding plane in the tissuespace, then the 2-D model is considered to be precise enough [1].

[0057] Thus, one can see that the RF image can be considered to be adistorted version of the true tissue reflectivity image, where thedistorting kernel is the transducer PSF and, consequently, the imagequality fairly depends on its spatial-frequency properties. The lateralsystem resolution (lateral “smearing” of the PSF) is inverselyproportional to the beam width, mainly depends on the aperture geometry.In order to obtain a uniform focus width throughout the image field,expanding apertures are used, in which the number of the receivingelements is reduced for shallow depth and gradually increased withdepth. Thus, one can assume the lateral resolution to be uniform alongthe image depth. The axial transducer resolution is a function of thewavelength and “time smearing” of PSF. The quality of ultrasound systemsis fairly defined by these resolutions.

[0058] This paper presents an algorithm for improving the axialresolution of the imaging system by the deconvolution of the RF image inthe axial direction. In such a case, we further simplify theabove-mentioned convolution equation and model the RF-image as combinedof 1-D RF-sequences. Such approach has a few drawbacks, one of which isthat the high correlation of the adjacent RF-sequences is not treated,appearing as a lateral smearing of the image. Yet for some applications,the proposed modeling may be quite sufficient and fruitful.

[0059] In the given model each RF-sequence is assumed to be aconvolution product of the 1-D PSF of the ultrasound transducer, s[f],and the spatial tissue response taken along one direction, f[n]2[5].Further, we will refer to the 1-D transducer PSF as an ultrasonic pulse,meaning that if the piezoelectric element electromechanical response wasa digital delta function, the 1-D transducer PSF might be viewed as somehypothetical finite length pulse sequence. The spatial impulse responsetaken along one direction will be referred to as a one dimensionalreflectivity function. Assuming scatterers on each RF-sequence arelocated on a uniform grid and the system impulse response (ultrasonicpulse) is shift invariant along each RF-line, the beam former gives asingle receive signal that is given by

g[n]=s[n]*f[n]+u[n]  (1.2)

[0060] Here u[n] is the additive white noise perturbing the received RFsignal g[n]. The influence of the convolutional kernel s[n] results in“deformation” of the useful information about f[n] and it is one of thefactors giving rise to appearance of the speckle noise. The specklenoise affects the quality ultrasonic images in which it is present byreducing the contrast, lowering the signal to noise ratio (SNR) andobscuring important diagnostic details. Thus, the task of recovering ofthe reflectivity function f[n] from the received signal g[n] is of greatimportance. Generally speaking, this problem is more complicate, becausethe system impulse response s[n] cannot be considered to be available.The absorption of the ultrasound energy in tissues increases with thepulse frequency. This frequency-dependent attenuation process causes toboth the pulse amplitude and shape to be dependent on depth in tissueand leads to the observed non-stationarity of RF-sequences. So, theproblem of the reflectivity function derivation or, the same, theproblem of tissue characterization can be reformulated in term of theblind deconvolution problem. One way to deal with the non-stationarityis to break up the RF-sequence into a number of (possibly overlapping)segments, such that it will be possible to state that within each onethe frequency-dependent attenuation process does not affect the pulsespectrum properties considerably [3]. In such a case it is possible totreat a given RF-sequence separately by the performing the deconvolutionprocedure for each segment. Finally, the deconvolved segments can beintegrated together supplying us the entire processed RF-sequence. Aproblem that could arise in such approach is that short data intervalsmay not contain enough data for precise estimation. The pulse estimationmethod proposed in this paper is tolerative enough to the data pointinsufficiency problem, and so, it was decided to use the above-mentioned“segment-wise” processing.

[0061] Our method is based on the homomorphic analysis of theRF-sequence spectrum [7], but in the proposed algorithm we do not usethe wide-accepted cepstrum-based technique [2, 4-6]. It has been shown[7] that for any signal of finite duration its complex cepstrum alwaysextend to infinity. Consequently, using the discrete Fourier transformfor its calculation always gives rise to aliasing errors. The cepstrumaliasing can be reduced (but not avoided) by zero padding of theRF-sequences. It also can be noted, that the setting the high samplingrate (it is intended to decrease aliasing errors in the reflectivityfunction spectrum estimation), increases the degrading effect of themeasurement noise on the cepstrum estimates. The higher the samplingrate, the wider the spectrum band, occupied almost exclusively by thenoise. This band has to be included to the cepstrum computation and,consequently, it causes significant errors in the estimation of thecepstrum and unavoidably in the estimation of the ultrasonic pulse.Another inconvenience is the fact, that for signal lengths about 512sample points (that is usually used for such algorithms) the probabilitythat most complex zeroes will be located close to the complex unitcircle is large, and it leads to instability of all cepstrum-basedmethods. To remove this instability exponential weighting with a factorα<1 is used, but unfortunately such weighting can give artifacts in theestimation of the ultrasonic pulse function.

[0062] Computation of the 1-D ideal low-pass filter parameters in thecepstrum domain also can be considered to be a quite problematic task,because it has to be changed systematically from some initial estimateto find optimal values giving the best deconvolution results for a givenprobe.

[0063] All the above-discussed problems related to the cepstrum basedprocessing brought us to another method for the ultrasound pulseestimation. The ultrasound pulse spectrum is recovered using the Waveletanalysis of the RF-sequence log-spectrum and reconstructed from theWavelet transform coefficients. After the pulse estimation the problem,expressed by (1.2), is reduced to a deconvolution problem as a specialcase of the inverse problem. This inverse problem can be solved in theLS sense, but its direct solution is impossible due to the spectralproperties of the autocorrelation matrix to be inversed. Thus, someregularization technique should be used. In order to implement suchinverse, instead of estimating the f[n], we find its stable moment(approximation) <f,w>, where w is a suitable mollifier, that can be, inWavelet language, a scaling function or a wavelet. Such approach iscalled the approximate inverse and it is explicitly explained in [10].

[0064] The proposed blind deconvolution algorithm based on the Wavelettransform technique is tested in 1-D deconvolution of computer simulateddata and in vitro ultrasound data in the radial direction. It was shownthat its performance remains stable for severe noise levels and obtainedresults are very correlative with the tested reflectivity functions.

[0065] Section 2 gives a brief discussion of the signals and itsspectra. A surface description of the homomorphic analysis is also givenhere. Section 3 exclusively deals with the main ideas of the Waveletanalysis and its use for the non-parametric pulse estimation algorithmderivation. Section 4 gives an overview on ill-posed inverse problemsand introduces a method for construction of the approximate inversefilter.

[0066] Finally, Section 5 demonstrates results of using the approach torecovering computer simulated reflectivity functions and reflectivityfunctions of ultrasound phantoms.

[0067] 2. Ultrasound Signals and Homomorphic Analysis

[0068] The frequency-dependent attenuation process appears as decreasewith distance of the mean frequency and the amplitude of the ultrasonicpulse. So, it is unacceptable to estimate the pulse using the wholedata, as in this case it may lead to serious estimation errors. It mightbe appropriate to break up the RF-sequence into some number of (possiblyoverlapping) segments and to state that within each one thefrequency-dependent attenuation process does not affect the pulsespectrum properties strongly. In this light, we can implement thedeconvolution procedure inside of each segment. The final result (i.e.,processed image line) might be seen as an integration of the segmentsafter the treatment. The model of the entire RF-line [3] can be asfollows $\begin{matrix}{{{g\lbrack n\rbrack} = {{\sum\limits_{k = 1}^{M}{{f\left\lbrack {n - k} \right\rbrack}{s\left\lbrack {k,{n - k}} \right\rbrack}}} + {u\lbrack n\rbrack}}},{n = 1},2,\quad \ldots \quad,N} & (2.1)\end{matrix}$

[0069] In this expression g(n) is a received RF-sequence, f[n] is areflectivity function, u[n] is the additive white noise term, and s[k,n]is the ultrasonic pulse ascribed to the point n. It was assumed thatdespite of the pulse dependence on the RF-line segment “depth”, it couldbe modeled as M-points finite length sequence. The finality of the abovesum also assumes that f(n) is a sequence with length no longer than Npoints.

[0070] To turn the problem to be locally stationary, we assume that thepulse statistic is unchangeable on intervals (n_(i−1), n_(i)], for i=1,2, . . . , I, n₀=0, n₁=N . It says that s_(i)(k)=s(k,n),n_(i−1)<n<n_(i). We set${f_{i}\lbrack n\rbrack} = \left\{ {\begin{matrix}{{f\lbrack n\rbrack},{n_{i - 1} < n < n_{i}}} \\{0,{otherwise}}\end{matrix}.} \right.$

[0071] For these assumptions the model (2.1) can be written as$\begin{matrix}{{g\lbrack n\rbrack} = {{\sum\limits_{i = 1}^{I}\quad {\sum\limits_{k = 0}^{M}\quad {{f_{i}\left\lbrack {n - k} \right\rbrack}{s_{i}\lbrack k\rbrack}}}} + {{u\lbrack n\rbrack}.}}} & (2.2)\end{matrix}$

[0072] So, we start the algorithm derivation for one convolution segmentof the above sum that can be modeled by the equation (1.2). Thisequation can be specified in the frequency domain by taking the Fouriertransform on both its sides.

G(ω)=F(ω)S(ω)  (2.3)

[0073] There are certain interesting observations that can be made for(2.3). Of interest here is the fact that the RF-sequence spectrum isobtained as a result of the multiplying the pulse spectrum by thereflectivity function spectrum and they are quite different regardingits statistical properties. One of the most acceptable models for theultrasonic pulse sequence is such in which it is defined as a smoothGaussian-like function multiplied by an exponential of certain phase[11]. The frequency-dependent attenuation process influences the pulsespectrum, and consequently, its time pattern, but it is reasonable toexpect that this influence is smooth with depth. It implies that,regardless the reflection depth, the pulse spectrum should be a smooth,“slow” function with a maximum in vicinity of the nominal resonancefrequency of the transducer. Such a spectrum cannot admit considerableruptures or spikiness.

[0074] In FIG. 1 one can see the pulse-echo signal received from aplanar reflector in a water tank. Assuming that the reflectivityfunction of this elementary reflector is a zero sequence with only oneweighted digital delta with appearance at the point appropriate to the“two-flight” time, the echo signal is to be close to our ultrasonicpulse.

[0075] One can see that the results presented in FIG. 1 are in a goodagreement with the above mentioned pulse model and the conclusions aboutthe pulse spectrum character. In order to conduct these measurements theultrasonic transducer and quadrilateral thick perspex stick used as aplanar reflector were immersed into a water tank. The transducer was thestandard case style Panametrics V384 3.5/0.25 MHz with diameter of0.46″. The Panametrics computer controlled pulser-receiver (model 5800)was used to excite it. The target was located in the distance lightlylonger than the transducer nominal focus length (4.515″). The obtainedechoes were received and digitized at 50 MHz using the Tektronix TDS-420digital oscilloscope. The other half of the story is the reflectivityfunction and its spectrum. It is well understood that the reflectivityfunction ascribed to speckled regions (soft tissue) of ultrasound imagesmay be an irregular, spiky, and although a normal distributed randomfunction exhibiting complexity of the underlying tissue structure [5].Another approach to the tissue reflectivity modeling is to assume thatit may admit some form of layered structure. This implies that thereflectivity function will be sparse, i.e. only a limited number ofsamples have nonzero value [3]. Let us to omit the discussion of whatmodel is more appropriate, because of importance here is that there is acommon point for both. Such reflectivity sequences will have veryirregular, spiky spectrum with numerous ruptures and dips. Moreover sucha spectrum should be broadband in contrast to the pulse spectrum. Hencethe RF-sequence spectrum is nothing more than a result of themultiplication of the irregular, “wild”, broadband reflectivity spectrumby the relatively smooth and band-limited spectrum of the ultrasonicpulse. In FIG. 2 one can see a simulated spectrum of the reflectivityfunction and a typical spectrum of the measured RF-sequence taken inlogarithmic scale.

[0076] Now, let us consider the complex logarithm of the G(w) in (2.3).The Fourier transform of the RF-sequence (i.e., Z-transform evaluated onthe unit circle) can be defined in terms of its amplitude and argumentG(w)=|G(w)|·e^(j·arg(G(w))). In the same manner, let us defineS(w)=|S(w)|·e^(j·arg(S(w))) and F(w)=|F(w)|·e^(j·arg(F(w))). Then thecomplex logarithm of the G(w) is defined as

log{G(ω)}=log{|G(ω)|}+j arg[G(ω)]=log{S(ω)}+log{F(ω)}==log{|S(ω)|}+jarg[S(ω)]+log{|F(ω)|}+j arg[F(ω)]  (2.4)

[0077] It implies that

log{|G(ω)|}=log{|S(ω)|}+log{|F(ω)|}arg[G(ω)]=arg[S(ω)]+arg[F(ω)]  (2.5)

[0078] The above relationships (2.4) and (2.5) are well known in thehomomorphic signal processing. Now, one can see that the log-amplitudeof the G(w) can be represented as a sum of two components, namely, thesmooth and regular log spectrum of the ultrasonic pulse and very jaggedand irregular log-spectrum of the reflectivity function. Moreover, if weglance at these functions from a point of its scales, it will be clearthat log{|S(ω)|} is much “slower”, than log{|F(ω)|} and it implies thatits scales are quite different. Hence, it is reasonable to assume thatits frequency spectra (i.e., frequency spectra of the log-spectra) donot overlap significantly and one can separate the two components orperform independent processing of the two components. Moreover, a lowfrequency pole (slowly varying component) and a high frequency plateau(rapidly varying component) typically characterize the Fourier transformof log{|G(ω)|}. This fact is used in the cepstrum-based methods for thedeconvolution, whose comprehensive description can be found in [4, 7].It already was pointed out in Section 1, that there are a number ofinconveniences related to it. The main problems are intrinsic infinityof the cepstrum leading to aliasing errors when using discrete Fouriertransform for its computation; sensitivity of the cepstrum to theout-of-band noise; necessity to compute the minimum phase and maximumphase cutoff values for the ideal filter in the cepstrum domain. Allthis brought us to proposal of another methods based on the principlesof the homomorphic deconvolution without entering the cepstrum domain.

[0079] Thus, to estimate the pulse log-spectrum we must to filter outthe irregular jagged component of the RF-sequence log-spectrum and forthis purpose the Wavelet analysis was used. The wavelet method “acts asa mathematical microscope in which we can observe different scales ofthe signal by just adjusting the focus”. In the given case we arelooking for a long-term, slow portion of the log spectrum. Such changesmight be better analyzed with a smooth, regular wavelet. In this workthe “Coiflet” family of wavelets have been employed.

[0080] 3. Main Ideas of the Wavelet Analysis and its Use for the PulseEstimation

[0081] We are far from the introducing this theory, but let us justhighlight a few basic points of this powerful, beautifully structuredanalysis, in order to supply an intuitive platform for the derivation ofthe pulse estimate algorithm.

[0082] The wavelet transform is defined by decomposing the signal into afamily of functions, which are the translations, and dilations of aunique function ψ(x). The function ψ(x) is called a wavelet and thecorresponding wavelet family is given by {{square root}{square root over(s)}ψ(s(s−u))}, where (s,u)ε

² [12]. The wavelet transform of a function f(x)εL²(

) is defined by $\begin{matrix}{{W\quad {f\left( {s,u} \right)}} = {\int_{- \infty}^{\infty}{{f(x)}\sqrt{s}{\psi \left( {s\left( {x - u} \right)} \right)}\quad {x}}}} & (3.1)\end{matrix}$

[0083] In order to reconstruct the analyzed signal from its Wavelettransform, the Fourier transform of the wavelet function should satisfythe following admissibility condition $\begin{matrix}{C_{\psi} = {{\int_{0}^{\infty}{\frac{{{\hat{\psi}(\omega)}}^{2}}{\omega}{\omega}}} < \infty}} & (3.2)\end{matrix}$

[0084] It implies that {circumflex over (ψ)}(0)=0, and that {circumflexover (ψ)}(ω) is small enough in vicinity of ω=0. If we denote {tildeover (ψ)}_(s)(x)=ψ_(s)(−x), where ψ_(s)(x)={square root}{square rootover (s)}ψ(sx), then we can rewrite the Wavelet transform at a point uand a scale s as a convolution product with {tilde over (ψ)}_(s)(x)

Wf(s,u)=f*{tilde over (ψ)} _(s)(u)  (3.3)

[0085] The resolution of the Wavelet transform varies with the scaleparameter s. The Wavelet transform decomposes the signal into a set offrequency bands having a constant size on a logarithmic scale [13]. Thesmaller the scale, the courser the resolution in the spatial domain andfine in the frequency domain, and vice versa. We can derive that thereconstruction of f(x) from Wf(s,u) is given by $\begin{matrix}{{f(x)} = {\frac{I}{C_{\psi}}{\int_{- \infty}^{\infty}{\int_{0}^{\infty}{{{Wf}\left( {s,u} \right)}{\psi_{s}\left( {x - u} \right)}{s}{u}}}}}} & (3.4)\end{matrix}$

[0086] Like a window Fourier transform a Wavelet transform is redundant.In other words, its value at a given point of the phase (time-frequency)space depends on its values at neighboring points of the phase space.

[0087] In order to overpower the transform abundance, the Wavelettransform can be discretized by sampling both the scale parameter andthe translation parameter. By choosing the scale parameter sampling tobe exponential, and the translation parameter sampling to beproportional to the former, the discrete Wavelet transform is defined by$\begin{matrix}\begin{matrix}{{W_{d}{f\left( {j,n} \right)}} = {{{Wf}\left( {a^{j},\frac{n\quad b}{a^{j}}} \right)} = {{\int_{- \infty}^{\infty}{{f(x)}{\psi_{a^{j}}\left( {x - \frac{n\quad b}{a^{j}}} \right)}d\quad x}} =}}} \\{{= {f*{{\overset{\sim}{\psi}}_{a^{j}}\left( \frac{n\quad b}{a^{j}} \right)}}},\quad {a - {{delation}\quad {step}}}}\end{matrix} & (3.5)\end{matrix}$

[0088] The choosing the scaling parameter a to be 2 and b to be 1provides a very important particular case of the Discrete Wavelettransform. It was proved that there exist some wavelets ψ(x) such that({square root}{square root over (2^(j))}ψ(2^(j)(x−2^(−j)))), where(j,n)εZ², is an orthonormal basis of L² (

). These particular wavelets are called orthogonal wavelets. In suchcase, any function can be reconstructed from its decomposition into aWavelet orthonormal basis with the classical expansion formula of avector into an orthonormal basis $\begin{matrix}{{f(x)} = {\sum\limits_{j \in Z}{\sum\limits_{n \in Z}{{\langle{{f(u)},{\psi_{2^{j}}\left( {u - {n\quad 2^{- j}}} \right)}}\rangle}{{\psi_{2^{j}}\left( {u - {n\quad 2^{- j}}} \right)}.}}}}} & (3.6)\end{matrix}$

[0089] One can define the vector space V₂ _(^(j)) composed from the setof all possible approximations of functions at the resolution 2^(j). Thesequence of vector spaces (V₂ _(^(j)) )_(jεZ) is called amultiresolution approximation of L₂(

). For any function f(x) at the resolution 2^(j) is given by theorthogonal projection of f(x) on (V₂ _(^(j)) )_(jεZ). An orthonormalbasis of V₂ _(^(j)) can be built by dilating and translating aparticular function φ(x) called a scaling function. The family offunctions (φ_(d) _(^(j)) (x−2^(−j)n))_(nεZ) is then the orthonormalbasis of V₂ _(^(j)) . The function φ(x) can be viewed as a low passfilter; therefore the discrete approximation of f(x) (i.e. itsorthogonal projection into (V₂ _(^(j)) )_(jεZ)) can be computed as

A ₂ _(^(j)) f=(f*{tilde over (φ)} ₂ _(^(j)) (2^(−j) n))=(<f(x),φ₂_(^(j)) (x−2^(−j) n>)_(nεZ)  (3.7)

[0090] Here the function {tilde over (φ)}(x) is a reflection of thefunction φ(x).

[0091] The difference of information between the approximations as theresolutions 2^(j) and 2^(j+1) (i.e., the details of the function at theresolution 2^(j)) are equal to the orthogonal projection of f(x) on theorthogonal complement of V₂ _(^(j)) in V₂ _(^(j+1)) . Let us denote thiscomplement space as O₂ _(^(j)) . The family of functions (ψ₂ _(^(j))(x−2^(−j)n))_(nεZ) is the orthonormal basis of O₂ _(^(j)) . When theresolution 2^(j) varies in the interval (0, ∞) the functions (ψ₂ _(^(j))(x−2^(−j)n))_(nεZ) constitutes a wavelet orthonormal basis of L₂(

) [12]. The difference of information between the approximations as theresolutions 2^(j) and 2^(j+1) can be computed by expanding the f(x) intothe orthonormal basis of O₂ _(^(j)) and it is characterized by the setof inner products

D ₂ _(^(j)) f=(f*{tilde over (ψ)} ₂ _(^(j)) (2^(−j) n))=(<f(x),ψ₂_(^(j)) (x−2^(−j) n>)_(nεZ)  (3.8)

[0092] The Fourier transforms of φ(x) and ψ(x) are characterized by thefollowing relations $\begin{matrix}\begin{matrix}{{{\hat{\varphi}(\omega)} = {\prod\limits_{k = 1}^{\infty}\quad {H\left( ^{{- }\quad 2^{- k}\omega} \right)}}},{{\hat{\psi}\left( {2\quad \omega} \right)} = {{G\left( ^{{- }\quad \omega} \right)}{\hat{\varphi}(\omega)}}}} \\{{{where}\quad {G\left( ^{{- }\quad \omega} \right)}} = {^{{- }\quad \omega}{H\left( ^{{- }\quad \omega} \right)}}}\end{matrix} & (3.9)\end{matrix}$

[0093] In the above expressions the filters G and H make a pair ofquadrature mirror filters. So, it is was proved, that A₂ _(^(i)) f andD₂ _(^(j)) f can be computed by a algorithm that is in fact a classicalscheme known in the signal processing community as the quadrature mirrorfilter bank decomposition or subband coding [14].

[0094] Now we can come back to the subject of our discussion. It wasshown that the log-spectrum of the RF-sequence is combined from theregular, smooth log-spectrum of the pulse and the irregular, jaggedlog-spectrum of the reflectivity function. Thus, it was concluded thatfrom the scale point of view these signals are quite different. As aconsequence of this fact it is reasonable to expect that there is aresolution level 2^(j) in which the approximation of the analyzed signalA₂ _(^(i)) f, where f=log{|G(w)|}, characterizes the low-frequencycomponents of f, dominated by log{|S(w)|}. In the same time, the detailsbear information, related to the jagged, high frequency components ofthe signal, induced by log{|F(w)|}. If we decide J to be the desiredresolution level, then the analyzed log-spectrum of the RF-sequence canbe decomposed as follows $\begin{matrix}{f = {{\log \left\{ {{G(w)}} \right\}} = {A_{2^{J}}f{\underset{j \leq J}{+ \sum}{D_{2^{j}}{f.}}}}}} & (3.10)\end{matrix}$

[0095] In the above expression we assume the pulse log-spectrumdominates in the first term, while the reflectivity log-spectrumdominates in the second.

[0096] In order to verify our assumption a 2⁸-point RF-sequence wassimulated according to (1.2). The pulse was defined as an exponentialwith a linear time-weighted Gaussian-like envelope and the reflectivitywas simulated as a function consisting of the digital delta functionswith random appearances and weights. Both the weights and theappearances were taken from the normal Gauss distribution. In FIG. 3 onecan see the generated function and the logarithm of its Fouriertransform. In order to decompose log{|G(ω)|} a wavelet filter of order 5from the Coiflet Wavelets family was chosen. It is a wavelet withsupport width of 29 points and the filter length of 30 points. Thiswavelet is an orthogonal compactly supported and it has a highest numberof vanishing moments for both the scaling and wavelet functions (i.e.,the vanishing moment's number for φ(x) is 30 and the vanishing moment'snumber for ψ(x) is 29). The maximal decomposition level was chosen to be4 and the decomposition was performed using the Fast Wavelet Transformalgorithm [14].

[0097] In FIG. 4 one can see results of the decomposition. Note, that atthis stage we intentionally did not add any noise to the simulatedsignal in order to clarify the results and, consequently, the main ideaof the approach. The noise will be added at the next stage and anappropriate method for its cancellation will be proposed.

[0098]FIG. 4 gives a clear picture of what happens in process of thedecomposition. It is obvious that the fourth level approximation of theRF-sequence log spectrum is quite smooth and regular. All the highfrequency information (i.e., ruptures, dips and spikiness of thespectrum) is included into the details of this decomposition. In orderto verify whether the obtained approximation is close to the logspectrum of the pulse, in FIG. 5 both spectra (i.e., estimated andoriginal) are plotted together. No doubts that our assumption is preciseenough. Thus, an important conclusion to be made here is that theamplitude of the ultrasonic pulse spectrum can be recovered from theapproximation coefficients of the Wavelet transform of the RF-sequencelog spectrum.

[0099] Now, we have to confirm the algorithm performance when the noisecomponent of the (2.1) is not zeroed. The simplest method to reduceout-of-band noise influence on the estimation process is to use aone-dimensional low-pass filtering of all RF sequences in the radialdirection. It is also possible to use the classical Wiener filter. InFIG. 6 the block-diagram of our process is viewed. We denote by g₀[n]RF-sequence before the noise addition.

[0100] It is easy to show that the frequency response H(ω) of theoptimal filter, with impulse response h(n), minimizing the varianceE{(g[n]−g₀[n])²} is given by [15] $\begin{matrix}{{H(\omega)} = \frac{{{S(\omega)}}^{2}}{{{S(\omega)}}^{2} + {\sigma_{u}^{2}/\sigma_{f}^{2}}}} & (3.11)\end{matrix}$

[0101] In the above expression S(ω) is the frequency response of theultrasonic pulse, σ_(u) ² is the noise variance, and σ_(f) ² is thereflectivity function variance. Note that it is valid in speckledregions of the RF-image, where the reflectivity function elements can bedefined as independently identically distributed. Otherwise, in place ofσ_(f) ² we have to use |F(ω)|², that is the power spectral density off[n].

[0102]FIG. 6: The Block Diagram of the Studied Process.

[0103] Since in the beginning of the processing S(ω) is unknown one candivide the filtering procedure into the two stages: deriving the firstestimate of the pulse spectrum from RF sequence filtered by a low-passfilter (out-of-band noise cancellation) by using the first stage resultsfiltering the RF sequence by the optimal Wiener filter given in (3.11)and, consequently, to compute the second, more precise estimate of thepulse spectrum

[0104] The second scheme performance might give us a more robustfiltering, but at the expense of greater computations. In ourexperiments, the former filtering scheme has been chosen.

[0105] Now we are about to get the pulse estimation in the time domain,but before it, the pulse Fourier transform phase must be computed. Incertain cases it is reasonable to assume the pulse is a minimum phasesequence. For rational systems it is equivalently to saying that allzeros and poles are inside the unit circle [7]. This requirement impliesthat the terms log|S(e^(jω))| and arg[S(e^(jω))] constitute the Hilberttransform pair, and as a consequence of the fact, arg[S(e^(jω))] can becomputed from log|S(e^(jω))|. Namely, knowing log|S(e^(jω))|, that isthe real part of the Fourier transform of the complex cepstrum of thepulse, allows us to compute the even (symmetric) part of the latter.This even part is the real cepstrum. If the ultrasonic pulse is minimumphase sequence, then its complex cepstrum decidedly real, causalsequence. Consequently, it completely defined by its even part, in otherwords, by its real cepstrum. So, in the case of minimum phase (ormaximum phase) pulses we can obtain the complex cepstrum from the realcepstrum, and practically it is computed by multiplying the realcepstrum with the following window $\begin{matrix}{{u_{+}\lbrack n\rbrack} = \left\{ \begin{matrix}{0,{n < 0}} \\{2,{n > 0}} \\{1,{n = 0}}\end{matrix} \right.} & (3.12)\end{matrix}$

[0106] At this point the reader could be a bit confused, because weclearly contradict ourselves. One of the motivations of the algorithmproposal was the desire to avoid possible aliasing errors in thecepstrum domain, and now we use it to reconstruct the pulse Fouriertransform phase. In our opinion, there is no significant contradictionhere, since log|S(e^(jω))| is quite smooth, “slow” function, and as aconsequence of the fact, its cepstrum is fairly concentrated in vicinityof the origin. Thus, it is unreasonable to expect here considerablealiasing errors.

[0107] Note that the assumption of the minimum phase pulse is very hardand might seem unrealistic. Yet it was proved that in many cases it isvalid and can be used for the pulse estimation purposes. Moreover, onecan use exponential weighting of the RF-sequence, and in such a case theconvolution in (1.2) is equal to a convolution of two sequences weightedexponentially. The zeros and poles of such a sequence are shiftedradially by the factor reciprocal to the exponential basis. Thus, theexponential weighting is a useful technique for converting a mixed phasesignal to a minimum-phase or a maximum-phase signal [7].

[0108] To verify the said above in FIG. 7 the ultrasonic pulse and itsminimum phase version are plotted together. It is obvious that theminimum phase assumption is not less basis. At this light, the proposedpulse estimation algorithm can be combined from two stages. At the firststage, we compute an approximation of amplitude of RF-sequence Fouriertransform and it supplies us an estimate of the pulse Fourier transformamplitude. At the second stage, the pulse complex cepstrum is computedand, consequently, the pulse is estimate in the time domain.

[0109] Now, we test the algorithm performance on synthetic datagenerated according to (2.1). The ultrasound pulse was minimum phase andits parameters were chosen as close as possible to ones of the measuredultrasound pulse (FIG. 7). The reflectivity was simulated as a functionconsisting of the digital delta functions with random appearances andweights. The signal-to-noise ratio was chosen to be approximately 7 dBand it is given by $\begin{matrix}{{SNR} = {10 \cdot {\log \left\lbrack \frac{\sigma_{signal}^{2}}{\sigma_{noise}^{2}} \right\rbrack}}} & (3.13)\end{matrix}$

[0110] Here σ_(signal) ² and σ_(noise) ² the signal and noise variances,respectively.

[0111] The estimation results are viewed in FIG. 8. On the left-handside picture one can see eight pulse estimates performed on eightdifferent data segments. On the right-hand side picture the originalpulse is viewed. Note that the estimation algorithm does not require anyinitial conditions or parameters, except of the analyzing wavelet andthe decomposition level choice. In FIG. 9 the original pulse is plottedvs. the averaged (from the previous eight pulses) estimated pulse.

[0112] The pulse estimation algorithm performance was tested on the realdata set. The RF-sequences were obtained using the standard case styletransducer (Panametrics V384 3.5/0.25 MHz with diameter of 0.46″). Theobtained echoes from a phantom, which specific structure will bedescribed in Chapter 5 of the paper, were received and digitized at 50MHz. An entire RF-sequence was divided into shorter sections. The numberof samples per section was 256. In FIG. 10 one can see ultrasound pulseestimated from five different RF-sections (left-hand side).

[0113] It can be seen that all the estimated pulses possess a similarform, resembling the ultrasound pulse, measured experimentally. Thevariations among the estimated pulses may be ascribed to the statisticalestimation errors, “coloration” of the reflectivity function (i.e., itmight not be statistically white sequence), and contribution from someeffects of the medium having spatially varying characteristics (such asaberration, for example). On the right-hand side of FIG. 10 averagedestimated pulse is viewed versus the measured ultrasound pulse, and itcan be seen that they are very close one to another. Additionally it isreasonable to assume that the pulse in adjacent lines is same, so anaverage pulse can be found for a number of lines. This averaging can beaccomplished either in the wavelet, Fourier or time domains and it isintended to decrease the estimate variance.

[0114] After the pulse estimation, the problem presented in (2.1), asthe blind deconvolution problem, becomes a classical inverse(deconvolution) problem. This procedure ultimately relies on knowledgeof the pulse and its preciseness. The pulse estimation stage wasperformed with moderate errors, and, consequently, it is reasonable toexcept that the following stage (deconvolution procedure) will be notdisturbed significantly by errors caused by the low acuteness of thefirst stage.

[0115] The next section deals exclusively with the inverse problemsolution. It will be shown how it is possible to stabilize its solutionby the approximate inverse technique. Results of simulations and invitro experiments will be shown also.

[0116] 4. Approximate Inverse for Ill-posed Problems and its Use for theReflectivity Function Estimation

[0117] Ignoring the noise term the equation (2.1) can be viewed as anoperator equation S f=g, where S is assumed to be linear operatorbetween the Hilbert spaces X and Y. In our case the operator S is theconvolution operator and its structure will be described below.

[0118] It is well understood that the multiplication of f by S is asmoothing operation that tends to dampen high-frequency components in f,such that the function g is a smoother function than f[16]. Thus, thedeconvolution process can be expected to amplify the high-frequencycomponents, and if the data was perturbed by noise, this deconvolutionprocess becomes highly unstable. The larger the frequency, the largerperturbations of the desired solution. So, regularizing techniques areneeded to stabilize the finding the optimal solution and it is a matterof this section. Singular value decomposition of square matrix A isdefined as $\begin{matrix}{S = {{U\quad \Lambda \quad V} = {\sum\limits_{i = 1}^{n}\quad {u_{i}\lambda_{i}v_{i}}}}} & (4.1)\end{matrix}$

[0119] The two matrices U and V are orthogonal (U^(T)U=V^(T)V=I) andconsist the left and right singular vectors. The middle matrix is adiagonal matrix and its diagonal elements are nonnegative singularvalues of matrix A ordered in non-increasing order. The well-knowndefinition of the condition number of S is given as $\begin{matrix}{{{S}_{2}{S^{- 1}}_{2}} = \frac{\lambda_{\max}}{\lambda_{\min}}} & (4.2)\end{matrix}$

[0120] Here λ_(max) and λ_(min) are maximal and minimal singular valuesrespectively. The well-known solution of the linear equation expressedin terms of the SVD is obtained as follows $\begin{matrix}{{g = {\sum\limits_{i = 1}^{n}{{\langle\quad {u_{i}^{T}g}\rangle}u_{i}}}},{f = {\sum\limits_{i = 1}^{n}\quad {{\langle{v_{i}^{T}f}\rangle}v_{i}}}}} & (4.3)\end{matrix}$

[0121] Thus we obtain $\begin{matrix}{f^{*} = {\sum{\frac{u_{i}^{T}g}{\lambda_{i}}v_{i}}}} & (4.4)\end{matrix}$

[0122] The singular values of the matrix S decay gradually to zero,consequently the condition number is approximately reciprocal of themachine precision, i.e. almost infinite, thus matrix S is singular.

[0123] The most prominent regularization method for ill-posed problemsis Tikhonov-Philips regularization, where the error norm is penalizedwith a priori information on the solution. For the sake of simplicity weuse the L₂-norm of the solution and get

ε² ={|S f−g∥ ₂ ² +γ∥f∥ ₂ ²}→min  (4.5)

[0124] The minimum of this functional is computed as solution of

(S ^(T) S+γI)f _(γ) =S ^(T) g  (4.6)

[0125] Parameter γ is a regularization parameter and when it approacheszero we obtain the original, non-regularized solution.

[0126] One can express f_(γ) in terms of the SVD of matrix S as follows${f_{\gamma} = {\sum\limits_{i = 1}^{n}{{F_{\gamma}\left( \lambda_{i} \right)}{\lambda_{i}^{- 1}\left( {u_{i}^{T}g} \right)}v_{i}}}},{{{where}\quad {F_{\gamma}\left( \lambda_{i} \right)}} = \frac{\lambda_{i}^{2}}{\lambda_{i}^{2} + \gamma}}$

[0127] The factors 0≦F_(γ)(λ_(i))≦1 are called the Tikhonov filterfactors and they control the damping of the individual SVD components ofthe regularized solution [17].

[0128] One can show that this regularization method is a special case ofthe approximate inverse. Approximate inverse means a solution operator,which maps the data g to a stable approximation of the solution of theill-posed problem S f=g. This inversion operator (stable inverse filter)can be precomputed without using g. Thus, we compute instead of f itsstable moment <f,e_(γ)>, where e_(λ) is a suitable mollifier. Thismollifier may be a low-pass filter, or in wavelet language a scalingfunction or a wavelet. Consequently, the high-frequency components of fwhich are mostly affected by noise are reduced. The computation of<f,e_(γ)> is achieved by approximating e_(λ) in the range of theoperator S^(T) by the reconstruction kernel (inverse filter) v_(γ),where S^(T)·v_(γ)=e_(γ). In such a case the following relationship takesplace

<f,e _(γ) >≈<f,S ^(T) v _(γ) >=<S f,v _(γ) >=<g,v _(γ)>  (4.7)

[0129] Here we assume the equation S^(T)v_(γ)=e_(γ) to be solvable, butif it is not, then one can approximate v_(γ) by minimizing the defect∥S^(T)v_(λ)−e_(γ)∥. This minimization problem leads to solution of thefollowing linear equation system

SS^(T)v_(γ)=Se_(γ)  (4.8)

[0130] There are diverse approaches to solve (4.8). It can be done bydirect inverse, or iteratively [16], or in the spectral domain [15]. Ofimportance here the fact that only light regularization (relativelysmall parameter γ for the Tikhonov approach, for example) is needed tofind the optimal v_(γ).

[0131] Now, let's come back to the deconvolution. If we solve ourdeconvolution task by minimizing the defect ∥S f−g∥², then it leads tothe following equation system

S^(T)Sf=S^(T)g  (4.9)

[0132] In the left-hand side of (4.9) the matrix S^(T)S is theautocorrelation matrix of the pulse, and it is positive semi-definiteToeplitz structured matrix. Thus, the left-hand side can be viewed as aconvolution of the reflectivity with the pulse autocorrelation function.Let us denote Õ=S^(T)S. In the right-hand side we clearly recognizematch filtering of the RF-sequence by the pulse. Let us denote b=S^(T)g.So, the problem to be solved is Õf=b and for obvious reasons it is illposed. In order to achieve its stable solution we use theabove-described approximate inverse approach. Note that in this case weseek for the inverse operator (filter) for the symmetric autocorrelationof the pulse.

[0133] As a suitable mollifier a low-pass filter with frequency response“surrounding” the pass-band pulse spectrum was chosen. The impulseresponse of such filter is the classical damped sinc-function. Thus,after the deconvolution (i.e., convolution with the inverse filter) weexpect to obtain a low pass filtered version (moment) of thereflectivity function. In FIG. 10 one can see the autocorrelationfunction of an estimated pulse and an approximate inverse filterappropriate to it.

[0134] Note the length of the inverse filter that is significant. Thelength depends on the spectral properties of the convolution kernel andthe noise level. In the given case the kernel, i.e., ultrasound pulse,is narrow-band, consequently, in order to achieve an acutedeconvolution, longer inverse filter should be used. The increasing thestabilization parameter of the inverse problem can reduce the filterlength, but it is immediately expressed in weaker deconvolution.

[0135] In the next section these methods will be applied to RF-sequencesmeasured in vitro from different tissue-mimicking phantoms.

[0136] 5. Experimental Results

[0137] Starting the experimental part of the work let us reminder thefollowing assumptions on which based the proposed pulse estimationprocedure.

[0138] The ultrasound pulse is deterministic, possibly minimum phase

[0139] The tissue response can be modeled as a stationary, white,zero-mean sequence, possibly Gaussian

[0140] The additive experimental noise is taken from the Normaldistribution with zero mean and variance σ_(u) ², and it isstatistically independent of the tissue response

[0141] If we process a RF-line segment belonging to the speckled regionof the interrogated tissue image, the above assumptions should hold.Outside such a region it is reasonable that the first and the thirdassumptions still exist, but the second could be violated. In order toavoid the estimate errors caused by the reflectivity function“coloration”, it is possible to employ any segmentation procedureextracting from an entire RF-line those segments, which seems to beobtained from randomly distributed scatterers.

[0142] Recordings: 12 RF-lines were recorded from the tissue-mimickingphantom (Multipurpose Tissue/Cyst Ultrasound Phantom, 84-317, NuclearAssociates Prod.) using standard case type transducer (PanametricsV384-SU) with central frequency of 3.5/0.25 MHz and the nominal focuszone of 4.515 inch. The transducer and the phantom were positionedinside a water tank, in such a way, that the target was located invicinity of the focus. Each RF-sequence consisted of 512 sample pointsobtained from the same depth of, approximately, 4-cm. The sampling ratewas chosen to be 50 MHz and the data acquisition was performed throughthe Tektronix-420 digital oscilloscope connected to a PC by a GPIBinterface. Note that much lower rates are usually used for thesepurposes, but this choice is justified by a necessity to avoid possiblealiasing errors in evaluation of the reflectivity function spectrum.

[0143] It was already pointed out, that the estimate errors are mainlycaused by the non-whiteness of the reflectivity function.Distance-dependent attenuation of the RF-signal appears as a component,giving rise to the “coloration” of the reflectivity function spectrum.To eliminate its negative influence, the measured RF-sequences weredivided into 2 shorter sequences of 256 sampling points each one. Notethat this procedure is also necessary for providing thequasi-stationarity condition, assuming that inside of short enoughRF-line segment the frequency-dependent attenuation does not distort thepulse spectral features strongly. Thus, the input data was structured inthe form of two 256-by-12 matrices, which are appropriate to different(but adjacent) depths. The pulse estimation algorithm was applied toeach column of these matrices, and subsequently, the obtained estimatedpulses (12 for each data matrix) were averaged to obtain one pulse foreach interrogated depth. Note that the averaging procedure can beimplemented in the Wavelet domain or in the log-spectrum domain, and itis intended to decrease the estimate variance (assuming that the pulsein the adjacent lines is same). It has been observed that the estimatedpulses obtained for the two data matrices (and so the appropriateaveraged pulses) are not considerably different and it seems to beconcentrated about a same mean value. It can be true, since it isreasonable to assume that the frequency-dependent attenuation processcannot appear significantly when comparing quite short (256 samplepoints) adjacent segments of the data. Analyzing the spectra of thesesequences also shows that they constitute same statistical properties.For these reasons, we would not like to demonstrate the results of thepulse estimation separately for each depth, but it is of interest tocompare the obtained estimation results with the ultrasound pulsemeasured in a water tank, when a thick planar perspex stick was used asa target. In the left-hand side of FIG. 12 eight pulse estimates areviewed, whereas the right-hand side of this picture compares its averagewith the measured pulse.

[0144] It is obvious that there is a reasonable match between theaverage estimated and the measured pulses, but it is also obvious thatin comparison with the measured pulse the estimated one seems to be“smeared”, i.e. it has lower central frequency and wider time support.To verify this fact it is instructive to compare between the spectrumsof the RF-sequence from which we estimate the spectrum of the pulse andthe measured pulse spectrum.

[0145] In FIG. 13 one can see the spectrum of one of the measuredRF-line (recording from the tissue mimicking phantom) superimposed bythe spectrum of the measured pulse (underwater recording from a thickplanar perspex stick). Obviously, the central frequency of the former isbiased toward lower frequencies and its “bell” seems to be a bitnarrower in comparison with the latter. Thus, the “smearing” of thepulse estimate is reasonable and can be accounted for by thefrequency-dependent attenuation processes existing in the phantommedium. Therefore, it leads us to the conclusion that the measured pulsedoes not truly represent the situation inside the issue-mimickingphantom.

[0146] Note that before the pulse estimation procedure, the inputsequences were filtered using a low-pass filter in order to remove theout-of-band noise influence. In our specific case central frequencies ofthe interrogating pulse are quite remote from the lowest frequencies,thus another filtering is needed to cancel the noise distorting thebeginning of the spectrum. Two standard LS designed filters were usedfor this purpose. The Wavelet decomposition up to level 4 was performedusing standard Mailat Fast Wavelet transform algorithm with the waveletfilter from the “Coiflet” family (N=5). Note that another waveletsfilter, such as minimum-phase wavelets of Daubechies (N=6, 8), symmetricwavelets of Daubechies (N=6), some biorthogonal wavelets alsodemonstrated a good performance. The interrogating pulse was presumed tobe minimum phase and its phase reconstruction was performed as describedin Chapter 3.

[0147] Now given the estimated ultrasound pulse we can construct anappropriate inverse filter and check its deconvolution performance. Toconstruct the inverse filter we have to choose a suitable mollifier, asit was described in Chapter 4. This choice depends mainly on the pulsespectrum and in the given case the scaling function constructed by Meyer(available in the Matlab toolbox through the routine “wavefun”), waschosen for these purposes. In the Fourier domain it constitutes alow-pass filter with a smooth, continuous cut-off slope [12], thatallows the effective support of the function in the time domain to bequite local (but generally speaking, the support of such a filter isinfinite). Note that using the ideal low-pass filter is never a goodchoice (despite it might seem reasonable from the spectral point ofview), since its impulse response extends to the infinity. Thus, bychoosing a differentiable spectrum (i.e., allowing the continuity in itscut-off region) it is possible to decrease significantly the negativeeffects of infinite filter lengths. Another choice of the mollifier isalso possible, and it is fairly a matter of opinion [10].

[0148] An additional difficulty encountered in the process of theinverse filter construction is caused by the ultrasound pulse spectralproperties. The used pulse was quite narrow band, and consequently theinverse filter is not able to “whiten” its spectrum perfectly (withinthe limits dictated by the mollifier pass-band). This situation isdepicted in FIG. 14, on the left-hand side of which one can see theinverse filter computed for the estimated pulse, and the right-hand sideshows the result of its convolution with the pulse itself. It should benoted here, that the inverse filter was computed for the pulse, in placeof its autocorrelation sequence, as it was described in Chapter 4, andonly insignificant correction is needed to do it. Regarding thedeconvolution results, we clearly see that the deconvolved resultresembles the typical low-pass filter impulse-response, but its sideripples are obviously larger in comparison with the latter. Thisartifact is fairly ascribed to the difficulty of the inversion ofnarrow-band spectra, and it may introduce an inaccuracy in the tissueresponse moment (low-frequency components) evaluation, and the latterappears as a weaker deconvolution. Yet, fortunately, in our specificcase the resulted kernel is close enough to be a low-pass filter impulseresponse, and no considerable errors are expected.

[0149] In despite of above-mentioned difficulties, it is now clear thatthe deconvolution of the RF-sequences should results in a significantreduction of the time-spatial support of the interrogating pulse,suggesting an apparent gain in the axial resolution. Therefore, we hopethat the deconvolved RF-lines will resemble low-pass filtered versionsof the appropriate tissue response sequences, and it will supply us morevalid information about the scatterers' location and distribution.

[0150] In order to investigate to performance of our technique, weapplied the deconvolution procedure to the RF-sequences obtained formthe same earlier used tissue-mimicking phantom. Starting the analysis ofthe obtained results we have to emphasize, that due to the finite sizeof the beam the received RF-sequences are not combined from echoes fromthe scatterers and reflectors located exclusively on a line that is anaxial symmetric axis of the transducer. A RF-line gets benefits from thescatterers and reflectors located in the entire 3-D region of theinsonified tissue defined by the beam thickness in the focus. Since ouralgorithm does not treat the lateral resolution, after the deconvolutionprocedure we cannot identify form which regions exactly (within thelimits of the ultrasound beam) the reflections were obtained. This factdisturbs us to check validity of the algorithm performance inexperiments with the tissue-mimicking phantom.

[0151]FIG. 15 (upper picture) shows the RF-sequence received from thetissue-mimicking phantom and filtered by an appropriate band-passfilter, whereas the lower picture demonstrates the same RF-line afterapplying the deconvolution procedure to it. Note that the insonifiedregion had no strong reflectors, therefore, it is reasonable that thisRF-line might contribute to a speckled region of an image. At the firstglance, the differences between the sequences are not so obvious, but ananalysis of its spectra shows that the deconvolved RF-line spectrum hasbeen whitened, i.e., it is almost flat in the frequency interval up tothe cut-off frequency of the noise canceling filter. Thus it leads us tothe conclusion that the obtained signal has the same statisticalproperties, which were assumed for the reflectivity function in thespeckle region. Unfortunately, it is not possible to prove whether thedeconvolved sequence indeed constitutes the low-pass filtered version ofthe original tissue response.

[0152] The idea to prove the validity of the scatterers' identificationbrought us to usage of a phantom specially constructed for this purpose.There are two principal features, which must be provided in constructionof such a phantom. First, it should be combined of parallel layersjointed together having a lateral size much larger in comparison withthe beam thickness. In such a case, we can ignore the artifactscontributed by non-treatment of the lateral resolution, and one can saythat the reflectivity function is now reasonably one dimensional, andall its scatterers are arranged along this dimension. Second, thedistances between the layers and its widths should be random (possiblytaken from the normal distribution). It also should be chosen to besmaller in comparison to the wavelength in order to supply us RF-lineswith significant echo interference.

[0153] This phantom is viewed in FIG. 16. It was constructed formseveral plastic and perspex sheets having different widths and jointedtogether in such a way that the distances between the sheets are randomwith a good precision. Despite of our attempts, for obvious reasons wecould not find plastics with depths taken from the normal distribution,and thus, only plastics with Δ=0.44, 0.5, 0.8, 0.9, 1.2, 1.44, 1.5, 1.8mm were randomly mixed and attached one to the another. Clearly, thereis one serious problem with such a phantom structure. Since the phantomis intended for underwater experiments, the pulse propagating through K,always passes sequentially from a more dense to a less dense, andsubsequently, from a less dense to a more dense matter (i.e., from theplastic to the water and vise versa). It other words the odd samples ofthe reflectivity sequence may be, for example, positive, whereas all theeven samples are negative. It implies that these samples are highlycorrelative, and it appears as a “coloration” of the spectrum. It wasalready pointed out that such deviation from the whiteness couldintroduce inaccuracy in the pulse estimation. Yet, fortunately, inprocess of the experimentation we observed that errors caused by thisproblem are not crucial.

[0154] Clearly that this phantom is far to be physical, but let'sreminder that it was created especially to check an axial resolutiongain after deconvolution. The pulse-estimating algorithm was applied tothe RF-sequences received from the phantom. It is reasonable to expectthat overall character of the pulse inside this specific phantom shouldnot be too unlike the pulse measured from a planar reflector orthogonalto the beam in the water tank. Thus, the capability of the algorithm toestimate pulses will be investigated by comparison of the estimatedpulses with results of the measurements. The results of the estimationalready have been shown in FIG. 10 of Chapter 3. In the right-hand sideof the picture the average estimated pulse is compared with the pulsemeasured in a water tank. Correspondence between the estimated and themeasured pulses is striking.

[0155] The estimated pulse was used to compute the inverse filter. Thisprocedure is absolutely analogous to one described earlier, and theinverse filter obtained in this case is negligibly different from thefilter viewed in FIG. 14 (right-hand side). The results of thedeconvolution are depicted in FIG. 17, where the upper picturedemonstrates a fragment of the measured RF-line, whereas the lowerpicture demonstrates a sequence obtained after applying the inversefilter to the former. The dashed line on both pictures is an envelope ofthe original RF-line, and it allows the visual comparison of thedeconvolved sequence with the conventionally processed sequence (dashedline). Visual comparison of the two plots makes additional explanationsredundant. Obviously, after the deconvolution we can easily distinctbetween the echoes reflected from surfaces of the phantom layers.Comparing it with the RF-envelope (dashed line), we conclude thatsignificant axial resolution enhancement is possible with deconvolution.

[0156] 6. Conclusions

[0157] This paper has demonstrated a novel pulse estimation techniquebased on the homomorphic signal processing and the Wavelet analysis ofthe log spectra. The proposed method is close in concept to thewell-known cepstrum-based methods [1], [2], [4-6], but there are alsoseveral conceptual differences, which makes the algorithm quiteattractive. The procedure is implemented in the log Fourier domainallowing to avoid some artifacts related to the cepstrum computation(aliasing error, for example). It was observed that the algorithm isquite tolerative to experimental noises, and its performance is quickenough (possibility of using in the real time) due to efficient scheme,known as Fast Wavelet transform, discovered by Mallat [14].

[0158] It was emphasized that not all the wavelet filters give identicalresults, and best performance is achieved with long, regular wavelets.This technique was applied to in vivo experimental data obtained fromthe different tissue mimicking phantoms, and it was shown that our pulseestimation procedure could be done with reasonable accuracy. Theaccuracy of the estimation was verified by underwater measuring theultrasound pulse from a thick planar reflector orthogonal to the beam.There is a reduction in the variance of the estimate, when severalestimated pulses are averaged.

[0159] In order to reconstruct phase of the pulses, minimum phaseassumption was made. Generally speaking, the interrogating ultrasoundpulse might not be minimum phase, but it was experimentally shown thatthe pulses used in the given case satisfy this condition with goodprecision. Currently, the topic of generalization of the algorithm forthe non-minimum phase conditions is an object of our intensiveattention. The estimated pulse was used in order to construct theapproximate inverse filter employed for deconvolution purposes.Deconvolution results obtained with the tissue-mimicking phantom dataand, in most degree, with the specially constructed “layers-build”phantom data indicate considerable axial resolution improvement. It hasbeen reported that wide-band pulses are preferable for achieving betterdeconvolution results.

[0160] Due to the apparent gain in the axial resolution after thedeconvolution those small structures (early stage tumor, for example)that are hidden away from view should be revealed. In addition to it, itwas shown that a successively deconvolved sequence should constitute amoment (low-pass filtered version) of the “true” reflectivity sequence(tissue response), that is largely independent of the imaging system. Asa consequence of the fact, these sequences might be used for tissuecharacterization schemes, which ideally require the above-mentionedindependence. It can be added that appropriate analysis of the pulseshape estimated from different depths of the interrogating tissue alsocould carry a useful information about the frequency deponentattenuation processes.

[0161] In conclusion, this paper establishes the 1-D deconvolution as apowerful technique removing the image distortion in the axial direction.Despite of the fact that such a procedure does not treat the “smearing”of images in the lateral direction, there are diverse kinds ofapplications where it could be very useful.

[0162] 7. Acknowledgment

[0163] We are very grateful to Zvi Fridman for providing the tissuemimicking phantom and his valuable suggestions, to Dan Bodik forconstructing the “layers-built” phantom. We wish to thank the anonymousreviewers whose comments substantially improved the paper.

[0164] References

[0165] [1] Udantha R. Abeyratne, Athina P. Petropulu, Jonh M. Reid,“Higher order spectra based deconvolution of ultrasound images”, IEEETransactions on Ultrasonics, Ferroelectrics and Frequency Control, vol.42, no. 6, November 1995

[0166] [2] Torfinn Taxt, “Restoration of medical ultrasound images usingtwo-dimensional homomorphic deconvolution”, IEEE Transactions onUltrasonics, Ferroelectrics and Frequency Control, vol. 42, no. 4, July1995

[0167] [3] Kjetil F. Kaaresen and Erik Bolviken, “Blind deconvolution ofultrasonic traces accounting for pulse variance”, IEEE Transactions onUltrasonics, Ferroelectrics and Frequency Control, vol. 46, no. 3, May1999

[0168] [4] Torfinn Taxt, “Comparison of cepstrum-based methods forradial blind deconvolution of ultrasound images”, IEEE Transactions onUltrasonics, Ferroelectrics and Frequency Control, vol. 44, no. 3, May1997

[0169] [5] J. A. Jensen and S. Leeman, “Nonparametric estimation ofultrasound pulses”, IEEE Transactions on Biomedical Engineering, vol.41, pp. 929-936, 1994

[0170] [6] J. A. Jensen, “Deconvolution of ultrasound images”,Ultrasonic Imaging, vol. 14, pp. 1-15, 1992

[0171] [7] A. V. Oppenheim and R. W. Schafer, Discrete time signalprocessing, London: Prentice Hall, 1989

[0172] [8] D. G. Childers, D. P. Skinner and R. C. Kemmerait, “Thecepstrum: A guide to processing”, Processing IEEE, vol. 65, pp.1428-1443, 1977

[0173] [9] T. J. Ulrych, “Application of homomorphic deconvolution toseismology”, Geophysics, vol. 36, pp. 650-660, 1971

[0174] [10] A. K. Louis, “Approximate inverse for linear and somenonlinear problems”, available fromhttp:/www.num.uni-sb.de/gebiete/invers/e_approx.html

[0175] [11] “Mathematical model for the ultrasound pulse sequence”

[0176] [12] I. Daubechies, “Different Perspectives on Wavelets”,American Mathematical Society, Short Course, Jan. 11-12, 1993

[0177] [13] Stephane G. Mallat, “Multifrequency channel decomposition ofimages and Wavelet models”, IEEE Transactions on Acoustic, Speech andSignal Processing, vol. 37, no. 12, May 1989

[0178] [14] Stephane G. Mallat, “A theory for multiresolution signaldecomposition: The Wavelet representation”, IEEE Trans. Pattern. Anal.Mach. Intell. 11, 1989, pp. 674-693

[0179] [15] E. Sekko, G. Thomas, A. Boukrouche, “A deconvolutiontechnique using optimal Wiener filtering and regularization”, Signalprocessing, 72, pp. 23-32, 1999

[0180] [16] Per Christian Hansen, “Numerical aspects of Deconvolution”,available from http/www.imm.tdu.dk/-pch

ANNEX B Robust Ultrasound Pulse Spectrum Estimation Algorithm UsingOutlier Resistant De-noising Dan Adam and Oleg Michailovich

[0181] 1. Speckle Noise

[0182] Due to coherence of the back-scattered echo signals, imagesobtained from echo ultrasound imaging systems have extremely complexpattern bearing no obvious relationship to the macroscopic properties ofthe insonified object. The vast majority of biological tissues areextremely rough on the scale of an acoustical wavelength, andconsequently a signal obtained within a resolution cell consists ofcontributions of many independent scatterers. Interference of thesede-phased echoes gives rise to pattern, which has appearance of achaotic jumble of “speckles”, known as speckle noise. The specklepattern consists of a multitude of bright spots where the interferencehas been highly constructive, dark spots where the interference has beenhighly destructive, and brightness levels between these extremes. Itappears chaotic and unordered, and it might be described quantitativelyby methods of probability and statistics. The presence of speckle in animage reduces the ability of a human observer to resolve fine details.It obscures very small structures (early stage tumor, for example)decreasing reliability of tissue characterization. Thus in most cases ofpractical interest, suppression of the speckle noise is a goal towardswhich we aspire.

[0183] 2. Basic Model

[0184] One approach of describing the scattering of tissue is toconsider the acoustic inhomogeneity as an assembly of reflectors andscatterers. The former is a plane interface large, compared to thewavelength, reflecting portions of the transmitted energy back towardthe transmitter. The latter is a small object, compared to thewavelength, spreading the transmitted signal in all directions,reflecting small portions of the energy. It is acceptable to model sucha system by a (generally speaking three-dimensional) function called thespatial response of insonified material, or in the case of medicalimaging spatial tissue response or reflectivity function. Below it willbe referred to as the reflectivity function.

[0185] One can consider an ultrasound radio frequency image (RF-image)as combined of 1-D echo-sequences, or using the terminology of theultrasound community, RF-lines. Assuming the tissue properties to beuniform in the plane perpendicular to the scanning beam, one can view anobtained 2-D RF-image as result of the convolution of the 2-Dreflectivity function (which accounts for inhomogeneity in the scanningplane) and 2-D transducer point spread function (PSF). Thus, theRF-image can be considered to be a distorted version of the truereflectivity function, where the distorting kernel is the transducerPSF. Its influence results in deforming the useful information about atissue (appearance of the speckle noise). The speckle noise affects thequality of ultrasonic images in which it is present by reducing thecontrast, lowering the signal to noise ratio (SNR) and obscuringimportant diagnostic details. So the task of recovering of thereflectivity function from the received signal is of great importance.Unfortunately, this problem is more complicate, since the system PSFcannot be considered to be available. The absorption of the ultrasoundenergy in tissues increases with frequency. This frequency-dependentattenuation process causes to both the pulse amplitude and shape to bedependent on depth in tissue and leads to the observed non-stationarityof RF-sequences.

[0186] 3. One-dimensional Approach

[0187] In medical ultrasound, a pulsed field is transmitted into thebody, and the echoes, back scattered toward the emitting transducer, aredetected as a voltage trace, RF-line. It has proved fruitful to modelthe RF-line directly as being result of combining a hypothetical 1-Dpulse with a hypothetical 1-D tissue. Assuming scatterers on each imageline be located on a uniform grid and the system impulse response isrange shift invariant along each image line, the beam former gives asingle receive signal equal to the coherent sum of the echoes only fromon grid scatterers. A discretized version of this receive signal can bewritten as $\begin{matrix}{{r\quad {f(t)}} = {{\sum\limits_{k = 1}^{N}{a_{k}{s\left( {t - k} \right)}}} + {n(t)}}} & (1.1)\end{matrix}$

[0188] Here rf(t) is RF-line, s(n) is the transmitted ultrasound pulse,k is twice the time-of-flight to the k-th reflector, and n(t) ismeasurement noise. Note that it is assumed here that the reflectorsthemselves do not modify the pulse except by a real scaling factora_(k). Using the convolution notation one can rewrite (1) as follows.

rf[n]=a[n]*s[n]+noise[n]  (1.2)

[0189] Here n may be referred to as a time index, and a[n] is areflectivity sequence. Since the frequency-dependent attenuation processappears as decrease with distance of the mean frequency and theamplitude of the pulse, it is commonly assumed that the receivedpulse-echo signal may be expressed as a depth dependent pulse convolutedwith tissue reflectivity function. Therefore, in (1.2) the pulse is tobe “location-dependent” (i.e. s[n] must be replaced by s[n, k], where kis the location index), and it leads to the observed non-stationarity ofthe RF-lines received from the tissue. In order to deal with thenon-stationarity it was proved to be fruitful to break up theRF-sequence into a number of possibly overlapping segments. In such acase it is reasonable to state that within each segment thefrequency-dependent attenuation process does not affect the pulsespectrum properties considerably and the equation (1.2) holds. So, theproblem of the reflectivity function derivation or, the same, theproblem of tissue characterization can be reformulated in term of theblind deconvolution problem, namely, for each segment of a given RF-linethe ultrasonic pulse should be estimated and “extracted” from it. Aproblem that could arise in such approach is that short data intervalsmight not contain enough data for precise estimation. The pulseestimation method proposed in this paper is tolerative enough to this“data insufficiency” problem, and thus, the above-mentioned“segment-wise” processing will be used throughout this paper.Additionally we note that a successively deconvolved RF-sequence shouldconstitute a moment (low-pass filtered version) of the “true”reflectivity sequence, that is largely independent of the imagingsystem. Consequently, these sequences might be used for tissuecharacterization schemes, which ideally require the above-mentionedindependence. In addition to it, due to the apparent gain in the axialresolution after the deconvolution those very small structures (earlystage tumor, for example), that are hidden away from view, should berevealed.

[0190] 4. Homomorphic Analysis and Spectral Separation

[0191] We start the algorithm derivation assuming that a receivedRF-line was segmented, and that a chosen segment to be deconvolved ismodeled by (2.1). This equation can be specified in the frequency domainby taking the Fourier transform on both its sides.

RF(w)=A(w)S(w)  (2.1)

[0192] Note that in (2.1) for the sake of simplicity we temporarilyignore the noise term. Of interest here is the fact that the RF-linespectrum is obtained as a result of multiplying the pulse spectrum bythe reflectivity function spectrum and they are quite different from theviewpoint of its spectral features. One of the most acceptable modelsfor the ultrasonic pulse sequence is such in which it is defined as asmooth Gaussian-like function modulated by an exponential of certainphase. It implies that the pulse spectrum should be a smooth, “slow”function with a maximum in vicinity of the nominal resonance frequencyof the transducer. Such a spectrum cannot admit considerable ruptures orspikiness. In FIG. 1 one can see the filtered pulse-echo signal receivedfrom a planar reflector in a water tank and its DFT (Discrete FourierTransform). Assuming that the reflectivity function of this elementaryreflector is a zero sequence with only one weighted digital deltaappearing at the point appropriate to the “two-flight” time, the echosignal is to be close to our pulse.

[0193] It is well understood that the reflectivity function ascribed tospeckled regions (soft tissue) of ultrasound images may be modeled as arandom process combined of independent identically distributed randomvariables exhibiting complexity of the underlying tissue structure. Suchreflectivity sequences will have very irregular, spiky DFT with numerousruptures and dips. Moreover such a spectrum should be broadband incontrast to the pulse spectrum. Hence the RF-sequence DFT is nothingmore than a result of the multiplication of the irregular, “wild”,broadband reflectivity function DFT by the smooth and band-limit DFT ofthe ultrasonic pulse. In FIG. 2 one can see a simulated periodogram ofthe reflectivity function and a typical spectrum of the measuredRF-sequence taken in logarithmic scale.

[0194] Now, let us consider the complex logarithm of the RF(w) in (2.1).The Fourier transform of the RF-sequence can be defined in terms of itsamplitude and argument RF(w)=|RF(w)|·e^(j·arg(RF(w))). Similarly, defineS(w)=|S(w)|·e^(j·arg(S(w))) and A(w)=|A(w)|·e^(j·arg(A(w))). As aresult, the following, well known in the homomorphic signal processingrelationships can be obtained.

log{|RF(w)|}=log{|S(w)|}+log{|A(w)|}arg[RF(w)]=arg[S(w)]+arg[A(w)]  (2.2)

[0195] We see that the log-amplitude of the RF(w) can be represented asa sum of two components, namely, the smooth and regular log-spectrum ofthe ultrasonic pulse and very jagged and irregular log-spectrum of thereflectivity function. If we glance at these functions from the point ofits scales, it will be clear that log{|S(w)|} is much “slower”, thanlog{|A(w)|} and it implies that its scales are quite different. Hence,it is reasonable to assume that its frequency spectra (i.e., frequencyspectra of the log-spectra) do not overlap significantly and one canseparate the two components. One approach for recovering thelog-spectrum of the pulse is based on the following fact. From themultiresolution analysis point of view, the pulse log-spectrum can beconsidered as a signal belonging to the subspace of the lowerresolution, while the reflectivity function log-spectrum is a signalwhich belong to the space of higher resolution. Therefore in order torecover the former all we need is to project the RF log-spectrum in anappropriate space. In order to implement the above-mentioned projection,the Wavelet transform is used. Adjusting the resolution level one canobtain an approximation of the analyzed signal, which in the given case,is a good estimate of the pulse log-spectrum. Clearly the pulse estimateitself can be obtained inverting the previous steps.

[0196] The above “projection” method is very fast and efficient.Technically it requires only wavelet decomposition of log{|RF(w)|}supplying us a vector of the wavelet coefficients, following by thelinear shrinkage (i.e., setting to zero) all the coefficients up to aprescribed level. Obviously, it is assumed that these shrunken (i.e.,zeroed) coefficients bear bulk of the energy associated with thelog-spectrum of the reflectivity function. After the linear shrinkage anestimate of the pulse log-spectrum is obtained by the inversetransformation. Despite its attractiveness this method demonstrates someinconvenience in optimal choices of the decomposing filter and thedecomposition level. Indeed, non-optimal choice of the level up to whichthe wavelet coefficients are to be shrunken may cause a sneaking ofwavelet coefficients induced by the “reflectivity” part of log{|RF(w)|}in the pulse log-spectrum estimate. Obviously, it might influencedestructively the estimate quality. On the other hand, if such a levelwere chosen “deep” enough, we would expunge some coefficients whosesignificance in the estimate is crucial. Moreover, the separability ofthe spectra in the Wavelet domain essentially depends on the filtersmoothness. Nevertheless, these algorithm parameters (i.e., thedecomposition “depth” and the filter regularity) can be optimized, andthe method demonstrates quite acceptable results. In the paper, anotherapproach for recovering the pulse log-spectrum will be proposed. Thismethod is non-parametric and much more robust. It considers log{|RF(w)|}as being combined from a signal to be recovered, i.e. the pulselog-spectrum and noise, which in the given case is the log-spectrum ofthe reflectivity function. Thus, in order to recover the log-spectrum ofthe pulse we have to “filter out” the irregular jagged componentsascribed to the latter. The procedure of recovering a desired signalfrom its measurements contaminated by noise using non-linear shrinkageof its wavelet coefficients is known as de-noising. In order toimplement it we will use the algorithm, which was proposed by David L.Donoho and Iain M. Johnstone in 1992. Note that in series of latterpapers, the authors significantly modernized their initial idea byproposing the adaptive de-noising (1993), ideal de-noising in a bestbasis chosen from a library of bases (1994), translation-invariantde-noising (in collaboration with Ronald R. Coifman, 1995). Thesealgorithms are consistent and progressive developments of the firstwork, and each work is a logical continuation of the previous one. Forthe reason of the space we don't test here performance of the pulseestimation using all these excellent de-noising methods, and only thefirst (i.e., initially proposed) algorithm will be used in the paper.

[0197] 5. Projection Method

[0198] We start by describing the earlier proposed method for recoveringthe pulse log-spectrum. In base of the algorithm is the property of theWavelet transform that causes the decomposition of a polynomial toproduce null details, when the number of vanishing moments of thewavelet exceeds the degree of the polynomial. Now we would like to citethe following theorem proved by I. Daubechies that is of greatimportance for our discussion.

[0199] Theorem 1:

[0200] If ψ_(m,n)(t)=2^(−m/2)ψ(2^(−m)t−n) constitute an orthonormal setin L²(

) with |ψ(t)|≦C·(1+|t|)^(−K−1−ε), ψ(t)εC^(K)(

) and ∂^(k)ψ/∂t^(k) is bounded for all k≦K, then

∫dt·t ^(k)·ψ(t)=0, for k=0, 1, . . . , K  (3.1)

[0201] From the above theorem stems that if analyzed signal ispolynomial with degree less that K, then its wavelet coefficients arezero for all n,m. It says that wavelets automatically suppress thepolynomials. Let us suppose that on the interval [a,b ] including 0, wehave expansion of the signal x(t) $\begin{matrix}{{x(t)} = {{x(0)} + {x \cdot \frac{\partial{x(0)}}{\partial t}} + {x^{2} \cdot \frac{\partial^{2}{x(0)}}{\partial t^{2}}} + \ldots + {x^{K} \cdot \frac{\partial^{K}{x(0)}}{\partial t^{K}}} + {y(t)}}} & (3.2)\end{matrix}$

[0202] In (3.2) the signal y(t) is the irregular part of the signalx(t). According the Theorem 1, the signals x(t) and y(t) have the samewavelet coefficients. In other words, the ψ wavelet suppresses theregular, polynomial part of x(t) and analyses the irregular part.

[0203] From the discussion of the previous section, the pulselog-spectrum is quite smooth signal, which can be considered as apolynomial of certain degree. Thus it is reasonable to expect that ifthe number of the vanishing moments of chosen wavelet exceeds thepolynomial degree, then its wavelet coefficients will be equal to zero.On the other hand, the reflectivity function log-spectrum is quiteirregular, noisy signal, which cannot be considered to be polynomial (atleast, polynomial with degree less than the vanishing moments number ofthe most wavelets in use), thus its wavelets coefficient are never zero.Consequently, the log-spectrum of RF-line can be considered combinedfrom the regular part (i.e., the pulse log-spectrum) and the irregularpart (i.e., the reflectivity function log-spectrum). Thus, if theanalyzing wavelet is chosen to have number of the vanishing momentsexceeding the polynomial order, then the wavelet coefficients are fairlyassociated with the reflectivity function log-spectrum, i.e. the waveletsuppresses the “slow” pulse log-spectrum which is locally assimilated toa polynomial.

[0204] Now let us describe the algorithm implementation from thepractical point of view. We define a family of DWT indexed by twointeger parameters L and M. Strictly speaking, regardingimplementational particularities we also have to define a kind of thetransformation, either “periodic” or “boundary adjusted”. For a fixedvalue of L and M we get a matrix W. Multiplication of data vector{y_(i)}_(i=0) ^(n−1), n=2^(J), by this matrix gives the waveletcoefficients of y.

w=W·y  (3.3)

[0205] It must be emphasized that in practice the DWT is not implementedby matrix multiplication, but by a sequence of finite-length filteringfollowed by a decimation step. The choice of the Wavelet transform isessentially a choice of filter. In the orthogonal case (which we choosefor simplicity in exposition) we have the inversion formula

y=W ^(T) ·w  (3.4)

[0206] Let w_(j,k) denote the wavelet coefficients indexed dyadically,such that j=0, . . . , J−1, k=0, . . . , 2^(j)−1 and W_(j,k) denote(j,k)-th row of W. Then the inversion formula (3.4) can be rewritten as$\begin{matrix}{y = {\sum\limits_{j,k}{w_{j,k} \cdot {W_{j,k}(i)}}}} & (3.5)\end{matrix}$

[0207] One can see that the data vector y can be represented as theabove linear combination of the basis element W_(j,k) with coefficientsw_(j,k).

[0208] For L<j<<n, 0<<k<<2^(j) we have the approximation

{square root}{square root over (n)}·W _(j,k)(i)≈2^(j/2)·ψ(2^(j) t),t=i/n−k·2^(−j)  (3.6)

[0209] Here ψ is the mother wavelet arising in a Wavelet transformationon

, as described in Daubechies (1988, 1992). She has shown that theparameter M controls the smoothness (number of the vanishing moments) ofψ. The greater the M, the smoother the mother wavelet.

[0210] Practically it is usual to perform the DWT (Discrete WaveletTransform) using the classical Fast Wavelet Transform (FWT) algorithm ofStephan G. Mallat using the filter-bank technique. This adapts to thepresent situation as follows. Let us denote an approximation of the datasignal at resolution level j=L by V_(L)y and details of the signal atresolution level j, or scale 2^(−j) by W_(j)y. It is given by$\begin{matrix}{{V_{L}y} = {{\sum\limits_{j < L}{{w_{j,k} \cdot W_{j,k}}\quad W_{j}y}} = {\sum\limits_{0 \leq {k\quad 2^{j}}}{w_{j,k} \cdot W_{j,k}}}}} & (3.7)\end{matrix}$

[0211] Then y can be recovered from these components via $\begin{matrix}{y = {{V_{L}y} + {\sum\limits_{L \leq j < n}{W_{j}y}}}} & (3.8)\end{matrix}$

[0212] Now suppose that a given RF-line log-spectrum was decomposedaccording to (3.8). If parameter M of the DWT was chosen sufficientlylarge, then the second term of (3.8) bears an information relatedexclusively to the reflectivity function log-spectrum. Thus,reconstruction from the only first term of (3.8), which is reasonablyconnected to the polynomial part of the signal under consideration andconstructed from the wavelet coefficients induced by the pulselog-spectrum will give us an acceptable estimate of the latter. It isclear, that the choice of the second decomposition parameter L is ofgreat importance, since “over-decomposition” might partially “expunge”coefficients induced by the desired signal to be recovered, and“under-decomposition” might cause leakage of undesired waveletcoefficients into the estimate. Yet, practically, there is no problem to“capture” an optimal L.

[0213] For the reason of the space we omit demonstration of performanceresults of the above-described “estimation-by-projection” algorithm,which may be found in [14]. And now we would like to introduce the novelalgorithm for the ultrasound pulse estimation based on the “de-noising”technique. Let us start this by a brief description of the de-noisingmethod proposed by David L. Donoho and Iain M. Johnstone (1992).

[0214] 6. De-noising by Soft-thresholding

[0215] In 1992 David L. Donoho and Iain M. Johnstone have proposed avery simple thresholding procedure for recovering functions from noisydata. It consists of three steps defined as follows.

[0216] First Step:

[0217] Data signal (preferably consisting of is n=2^(J) sample points,where n,JεZ) is transformed into an orthogonal domain using periodizedinterval-adapted pyramidal algorithm of Cohen, Daubechies, Jawerth andVial (1993), obtaining empirical wavelet coefficients.

[0218] Second Step:

[0219] The empirical wavelet coefficients are coordinatewise subjectedto the soft-thresholding non-linearity η₁(y)=sgh(y)·(|y|−t), withspecially chosen threshold t_(n) ={square root}{square root over(2·log(n))}·σ, where σ is the standard deviation of the white noise.

[0220] Third Step:

[0221] The thresholded wavelet coefficients are inversely transformedinto the original domain supplying us a noise-free estimation of thetrue signal. We are far from to repeat all mathematical rigors of theexperienced authors, which might be found in the original paper, butnevertheless, let us to point up several important features of thealgorithm demonstrating its powerfulness and uniqueness. We remind thatthe problem of the recovering of an unknown function from noisy data canbe formulated as follows. Let {d_(i)}_(i=0) ^(n−1) denotes the noisymeasurements of an unknown function f(t)on [0, 1]. Supposing the noiseto be normal, independent identically distributed we have the followingmodel $\begin{matrix}\begin{matrix}{{d_{i} = {{f\left( t_{i} \right)} + {\sigma \cdot z_{i}}}},{i = 0},\quad \ldots \quad,{n - 1}} \\{{{where}\quad t_{i}} = {{{i/n}\quad {and}\quad z_{i}^{i.i.d}} - {N\left( {0,1} \right)}}}\end{matrix} & (4.1)\end{matrix}$

[0222] The goal is to optimize the mean-squared error $\begin{matrix}{{{\frac{1}{n} \cdot E}{{\hat{f} - f}}_{1_{n}^{2}}^{2}} = {\frac{1}{n} \cdot {\sum\limits_{i = 0}^{n - 1}\quad {{E\left( {{\hat{f}\left( {i/n} \right)} - {f\left( {i/n} \right)}} \right)}^{2}.}}}} & (4.2)\end{matrix}$

[0223] It is important to emphasize that this minimization is subjectedto the following additional side condition

with high probability, {circumflex over (f)} is at least as smooth asf  (4.3)

[0224] In order to implement this constrained minimization task theabove-mentioned three-step algorithm was proposed. It was shown thatthis approach provides better quality than procedures based onmean-squared error alone.

[0225] It was completely proved that the proposed three-step algorithmhas two unprecedented properties. First, with high probability theestimated signal is at least as smooth as the original one, withsmoothness measured by any of a wide range of smoothness measures.Mathematically it is rephrased as the following. Let S to be the scaleof all spaces B_(p,g) ^(σ) (Besov classes) and all spaces F_(p,g) ^(σ)(Triebel classes) which embed continuously in C[0, 1], {{circumflex over(f)}_(i) ^(σ)(t_(i))}_(i=0) ^(n−1) be the vector of estimated functionvalues produced by the algorithm, and {circumflex over (f)}_(n) ^(σ)(t)is a special smooth interpolation of these values producing a functionon [0, 1]. There are universal constants (π_(n)) with π_(n)→1, n=2^(j)^(₁) →∞, and C₁(F,ψ), which depends on the function space F[0, 1]εS andon the wavelet basis, but not on n or f, so that

PROB{∥{circumflex over (f)} _(n) ^(σ)∥_(F) ≦C ₁ ·∥f∥ _(F) , ∀FεS}≧π_(n)  (4.4)

[0226] In words, {circumflex over (f)}_(n) ^(σ)(t) is, with overwhelmingprobability, simultaneously as smooth as f in every smoothness space Ftaken from the scale S.

[0227] The second amazing property of the proposed algorithm is that theestimate achieves almost the minimax error over every one of a widerange of smoothness classes. Formally, let F[0, 1] be a function space(one of the Besov or Triebel spaces) and let F_(C) denote the ball offunctions {f:∥f∥_(F)≦C}. Then the smallest error that can be achievedfor any estimator, uniformly over all fεF_(C) is the minimaxmean-squared error given by $\begin{matrix}{{M\left( {n,F_{C}} \right)} = {\inf\limits_{\hat{f}}\quad \sup\limits_{F_{C}}{\frac{1}{n} \cdot E}{{\hat{f_{n}^{*}} - f}}_{1_{n}^{2}}^{2}}} & (4.5)\end{matrix}$

[0228] For each ball F_(C) arising from FεS, there is a constantC₂(F_(C),ψ) which does not depend on n, such that for all n=2^(j) ^(₁),j₁>j₀, $\begin{matrix}{{\sup \quad \underset{f \in F_{C}}{E{{\hat{f_{n}^{*}} - f}}_{1_{n}^{2}}^{2}}} \leq {{C_{2} \cdot {\log (n)} \cdot \inf\limits_{\hat{f}}}\quad \sup\limits_{F_{C}}E{{\hat{f} - f}}_{1_{n}^{2}}^{2}}} & (4.6)\end{matrix}$

[0229] In words, {circumflex over (f)}_(n) ^(σ)(t) is simultaneouslywithin a logarithmic factor of minimax error over every Besov, Holder,Sobolev, and Triebel class contained in C[0, 1].

[0230] The properties are of great importance for understandingadvantages of the proposed algorithm. The first one says that arecovered function with high probability is not contaminated by noise,i.e. the algorithm rejects the noise completely. The second one saysthat the estimate is nearly minimax over very wide range of smoothnessclasses. In general, existing non-wavelet methods can achieve the samesuccess only over a limited range of the balls F_(C) arising in thescale S (basically L² Sobolev balls only), by relatively complicatedmeans

[0231] All above-mentioned advantages are really due to the Waveletbasis, and it was proved, that no analogous results could be achieved byusing thresholding (or other possibly non-linear techniques) in theFourier basis. The method works because the Wavelet transform of anoiseless (true) signal compresses its l² energy into a relative smallnumber of large coefficients. On the other hand, Gaussian white noise inany one orthogonal basis is again a white noise in any other and withthe same amplitude. Thus, in the orthogonal Wavelet basis, the “few”non-zero signal coefficients really stick up above the noise. Therefore,the thresholding of the wavelet coefficients has the effect that itkills the noise while not killing the signal. Let's define thecompression number c_(n) (it measures how well a vector f can beapproximated by a vector with n nonzero entries) as $\begin{matrix}{c_{n} \equiv {\sum\limits_{k > n}{f}_{(k)}^{2}}} & (4.7)\end{matrix}$

[0232] Comparing of the extents to which the energy is compressed into afew big coefficients in the Fourier and Wavelet domains it may bedemonstrated that in most practical cases for the same n the magnitudeof the c_(n) in the Wavelet domain greater or equal to that in theFourier domain. It guarantees that the ideal diagonal projector worksbetter in the Wavelet basis than in the Fourier basis. The compressionadvantages of the Wavelet basis are responsible for the mean-squarederror advantages of the thresholding of the wavelet coefficients.

[0233] Generally speaking the proposed de-noising algorithm is suited towork well, when a measured signal is contaminated by the white noise.Fortunately, it was showed that insignificant corrections are needed toadjust the algorithm to several other (not all) noise types. Returningto our estimation problem we remind that the measured signal is thelog-spectrum of a RF-line, and it is assumed combined from the desiredsignal to be estimated, i.e., the pulse log-spectrum, and noise, whichis here the reflectivity function log-spectrum. A question that shouldbe asked at this point is what's a kind of noise we have, and whetherthe de-noising algorithm is able to reject it? In order to answer thisquestion we have to discuss statistical properties of the reflectivityfunction log-spectrum.

[0234] 7. Reflectivity Log-spectrum Distribution

[0235] We assume that a chosen segment of RF-line obtained from a softtissue is a coherent sum of echoes reflected from greatly large numberof scatterers having randomly distributed reflection coefficients, andrandom, disordered locations. In such a case it is acceptable to thinkthat according to the large numbers low, the scatterers distribution isclose to be Gaussian, and consequently further we suppose that RF-lineobtained from soft tissues is convolution of the ultrasound pulse withthe white noise {a_(i)}_(i=1) ^(n)˜N(0,σ²). Let {ReA_(i)}_(i=1) ^(n)denote the real part of the DFT of {a_(i)}_(i=1) ^(n) and{ImA_(i)}_(i=1) ^(n) denote its imaginary part. Due to the orthogonalityproperties of the DFT {ReA_(i)}_(i=1) ^(n) and {ImA_(i)}_(i=1) ^(n) areuncorrelated and normally distributed. $\begin{matrix}{{{{Re}\quad A} = {{\frac{\left\{ {{Re}\quad A} \right\}_{i = 1}^{n}}{\sqrt{n/2}} \sim {{N\left( {0,\sigma^{2}} \right)}\quad {and}{\quad \quad}{Im}\quad A}} = {\frac{\left\{ {{Im}\quad A} \right\}_{i = 1}^{n}}{\sqrt{n/2}} \sim {N\left( {0,\sigma^{2}} \right)}}}}\quad} & (5.1)\end{matrix}$

[0236] It can be shown that the probability distribution of Z={squareroot}{square root over (ReA²+ImA²)} is the Rayleigh distribution whoseprobability density function is given as follows. $\begin{matrix}{{f_{z}(z)} = {{\frac{\quad}{z}{F(z)}} = \left\{ \begin{matrix}{{\frac{z}{\sigma^{2}}{\exp \left( {- \frac{z^{2}}{2\sigma^{2}}} \right)}},} & {z \geq 0} \\{0,} & {z < 0}\end{matrix} \right.}} & (5.2)\end{matrix}$

[0237] According to (2.2) we compute the logarithm of Z and we wish toinvestigate how this process affects its statistics. The densityfunction of Y=log(Z) is given by $\begin{matrix}{{{f_{Y}(y)} = \frac{f_{z}(z)}{{{y}/{z}}}}{{{Using}\quad \frac{y}{z}} = {\frac{1}{z} = {\frac{1}{\exp (y)}\quad {we}\quad {derive}}}}} & (5.3) \\{{f_{Y}(y)} = {{\frac{z^{2}}{\sigma^{2}} \cdot {\exp \left( {- \frac{z^{2}}{2\sigma^{2}}} \right)}} = {{2{\frac{\exp \left( {2\quad y} \right)}{2\sigma^{2}} \cdot {\exp \left( {- \frac{\exp \left( {2\quad y} \right)}{2\sigma^{2}}} \right)}}} = \quad {= {2 \cdot {\exp \left( {\left\lbrack {{2y} - {\ln \left( {2\sigma^{2}} \right)}} \right\rbrack - {\exp \left( \left\lbrack {{2y} - {\ln \left( {2\sigma^{2}} \right)}} \right\rbrack \right)}} \right)}}}}}} & (5.4)\end{matrix}$

[0238] The probability density function f_(γ)(y) is of the Fisher-Tippetform (or doubly exponential form). When the second exponential isreplaced by the first three terms of its serial expansion, we have$\begin{matrix}\begin{matrix}{{f_{Y}(y)} \cong {\frac{2}{e} \cdot {\exp \left( {- \frac{\left\lbrack {{2y} - {\ln \left( {2\sigma^{2}} \right)}} \right\rbrack^{2}}{2}} \right)}}} \\{= {\frac{2}{e} \cdot {\exp \left( {{- \frac{1}{2}} \cdot \frac{\left\lbrack {y - \frac{\ln \left( {2\sigma^{2}} \right)}{2}} \right\rbrack^{2}}{0.25}} \right)}}}\end{matrix} & (5.5)\end{matrix}$

[0239] The significance of this expression is that the reflectivityfunction log-spectrum approximately can be viewed as a gaussian, randomvariable with the mean {overscore (Y)}=1n(2σ²)/2 and the constantvariance σ_(γ) ²=0.5². It is quite amazing result saying that regardlessthe variance of the tissue reflectivity functions, its log-spectrum(i.e., logarithm of its DFT) is white noise with constant variance. Fromthe viewpoint of the de-noising, we have a problem to cancel the whitenoise whose variance never changes. Note that the variance of thereflectivity functions does influence the mean of the distribution. Yet,it is unnecessary to try to compute and subtract it from RF-linelog-spectrum, since its presence has no any destructive effect on theestimate quality (i.e., its only influence is multiplication of theestimate by a constant factor).

[0240] In FIG. 3 (upper plot) one can see histograms (continuous lines)of simulated reflectivity function log-spectra built for differentstandard deviations of the reflectivity sequences, namely 0.5, 5, 40,and 200 (from left to right in the picture). The dotted lines denote thetheoretically computed Fisher-Tippet probability density functions. Onecan see that, indeed, the mean values of the distributions depend on thereflectivity function variances ({overscore (Y)}=1n(2σ²)/2), whereas itsvariances remain approximately constant. FIG. 3 (lower plot) shows theFisher-Tippet distributions with its Gaussian approximations. One cansee that in vicinity of the mean values the approximations are quiteexact, whereas achieving the peripheries the approximation error becomesmore significant. Nevertheless, very nearly we can accept thisapproximation as quite satisfactory.

TABLE NO 1 Comparison of first two moments of the Fisher-Tippetdistribution and its Gaussian approximation computed for differentvalues of the tissue reflectivity sequences Reflectivity Fisher-TippetGaussian function distribution approximation standard deviation meanvariance mean variance 0.5 −0.62209 0.3736 −0.34657 0.25 5 1.66760.41002 1.956 0.25 40 3.7469 0.41 12 4.0355 0.25 200 5.3563 0.411235.6449 0.25

[0241] In the Table No1 the first two moments of the Fisher-Tippetdistribution and its Gaussian approximation are compared. We see thatindependently of the reflectivity function variance the variance of theformer is almost constant, that is in good agreement with the behaviorof the Gaussian approximation variance. TABLE NO 2 Comparison ofdifferences between first two moments of the Fisher-Tippet distributionand its Gaussian approximation computed for different values of thetissue reflectivity sequences Difference between the Difference betweenthe Reflectivity means of the Fisher-Tippet variances of the Fisher-function distribution and its Tippet distribution and standard Gaussianapproximations its Gaussian approximations deviation (absolute values)absolute relative (%) 0.5 0.2755 0.1236 33.08 5 0.2884 0.1600 39.03 400.2886 0.1612 39.20 200 0.2886 0.1612 39.21

[0242] Table No2 shows the absolute and relative deference between thefirst two moments of the Fisher-Tippet distribution and its Gaussianapproximation computed for different values of the reflectivity sequencevariance. We note that the difference is almost constant, i.e.independently of the reflectivity function variance the mean of theapproximation is shifted by a constant value relatively to the exactdistribution mean, and likewise the ratio between the approximationvariance and the real distribution variance is almost constant. From theabove numerical results it stems that roughly we can take the valueabout 0.25-40%, as a priori guess about the “noise” variance.

[0243] It should be emphasized that the above Gaussian approximation isonly an approximation, and frankly speaking there is a cardinaldifference between these distributions. Attention should be drawn to thefact that the Fisher-Tippet distribution is non-symmetric. Let ushighlight this point in more details.

[0244] In FIG. 4 the Fisher-Tippet distribution is shown at once withits Gaussian approximation. Its right-hand side has a form nearlyresembling the “Gaussian ball”, but obviously having slightly smallervariance. One can observe that in this part the Gaussian distributionapproximates quite satisfactory the exact one. In contrast the left-handside of the Fisher-Tippet has a quite long “tail”, giving rise to thevariance larger than 0.25. Obviously, our approximation falls righthere.

[0245]FIG. 5 shows a typical Gaussian white noise sequence with unitvariance (upper plot) and a typical Fisher-Tippet i.i.d. “noise” withvariance σ² equal to 0.25 (lower plot). Note the straight linesdetermining regions of two standard deviations for both distributions.One can see that majority of the white noise samples (realizations) isconcentrated in the region of ±294 and only its minor portion slightlyextends beyond this region. In the same time the Fisher-Tippet behavesdifferently. Majority of its realizations is located within ±2σ regionresembling the white noise behavior, but there are also a fewrealizations, which are significantly “remote” from it (the left-hand“tail” of the Fisher-Tippet distribution is “responsible” for itsappearance). Thus nearly this noise could be viewed as white noisecontaminated by occasional transients or “outliers”. It was mentionedearlier that the three-step de-noising algorithm (as it was formulatedinitially) is more suitable for rejection of white, gaussian noises,inasmuch as it is based on a fundamental statistic consequence of theorthogonality of the DWT that transforms white noise into white noise ofthe same level. Thus it is reasonable to assume that the “outliers”might play a dramatic role here. There are two questions we have to ask.First, may we still use the originally proposed formula for thesoft-threshold calculation, namely t_(n) ={square root}{square root over(2·log(n))}·σ (here n is a number of data points and σ is noise level)?Second, what is a better way to estimate the noise variance σ ²?

[0246] In the previous chapter we have discussed the probabilitydistribution of the “noise” to be rejected. Some conclusions were madeon a behavior of the first two moments of its distribution.Nevertheless, in practice we often prefer to estimate it from data. Notethat it does not contradict the above-obtained results, but we must notethat it was obtained regardless possible measurement noises, and inorder to turn the algorithm to be adaptive to possible “surprises”, itmake sense to estimate the moments. Usually, when a signal to bede-noised is band-limited and sampled by a frequency exceeding at leasttwice its Nyquist rate, the white noise variance can be estimated usingwavelet coefficients of the finest decomposition level. Indeed, in sucha case it is reasonable to assume that no induced-by-signal coefficientsare present in this level, and only induced-by-noise coefficients arecaptured here. Using the fact that the white noise remains a white noisewith the same level in the Wavelet domain, its variance can be estimatedfrom these empirical wavelet coefficients. However, in numberapplications it could not be guaranteed that coefficients of the finestlevel do not contain the signal coefficients at all. Such situation istypical for number of the image processing applications, where edgesusually induce wavelet coefficients in several levels including thefinest. In order to obtain unbiased noise variance estimate based on thenoise measurements in presence of few coefficients ascribed to theoriginal signal, the following robust estimate is used. We can estimatethe variance of the noise by assuming that “most” empirical waveletcoefficients in the finest of the decomposition are induced by thenoise, and thus the median absolute deviation reflects the size of thetypical noise. Dividing this value by 0.6745 we obtain a robust estimateof the noise standard deviation.

σ=MAD/0.6745, where MAD=median(|wc−median(wc)|).  (6.1)

[0247] Here n is the data length and wc is the wavelet coefficient inthe finest resolution levels.

[0248] Now let's return to our specific case and answer the questionabout a best strategy for the noise variance estimate. It has beenargued previously that the pulse log-spectrum can be locally assimilatedto a polynomial, thus it is reasonable to assume that the first (finest)decomposition level does not contain coefficients of the signal.Consequently, one can conclude that only coefficients we observe hereare of the noise. At first glance, such a situation does not requireusing the robust median estimate according to (6.1), and we can computeit classically as $\begin{matrix}{\sigma = \sqrt{\frac{1}{n - 1} \cdot {\sum\limits_{i = 1}^{n}\quad {{{wc} - {{mean}({wc})}}}^{2}}}} & (6.2)\end{matrix}$

[0249] Nevertheless, we prefer the median estimate for the followingreason. The structure of the i.i.d. Fisher-Tippet noise is such that itcan be approximately viewed as “contaminated-by-outliers white noise”(see above discussion about the Gaussian approximation of theFisher-Tippet distribution, non-symmetry of the latter, and see FIG. 5for illustration). It is clear that these “transients” are captured inthe relatively large wavelet coefficients of the finer resolutionlevels, whereas the other noise samples give rise to the waveletcoefficients distributed more or less with the same variance throughoutall the levels. In such a case an estimate computed according to (6.2)from coefficients of the finest level consisting of the coefficientsinduced by the “gaussian” part of the Fisher-Tippet noise and boundedwithin the region of ±2σ, and few large “induced-by-outliers”coefficients, will be significantly biased by the latter. This estimatefairly reflects behavior of the second moment of the noise, but it mightbe greatly “dangerous” in the following sense. The greater the variance,the larger the threshold according to t_(n) ={square root}{square rootover (2·log(n))}·σ. Surely, that such a threshold guarantees suppressionof the “induced-by-noise” coefficients, but also it might be too largeand, consequently, it might suppress somewhat “induced-by-signal”wavelet coefficients causing oversmoothing of the result. In otherwords, few extreme, “from-the-distribution tail” noise samples mightbias (enlarge) the variance estimate and, subsequently, the thresholdvalue significantly, causing distortion (oversmoothing) of the de-noisedsignal. To avoid such situation, it makes sense to use (6.1) ignoringthe “induced-by-transients” wavelet coefficients in the estimate. Ofcourse in such a case, the threshold will kill most the noisecoefficients, while allowing to the “induced-by-outliers” coefficientsto pass into the signal estimate. In other words, these coefficients areviewed as being ascribed to the desired signal and having significantenergy (due to the excellent compressing properties of the Wavelettransform) against a Gaussian white noise background. Subsequently, thethresholding procedure will keep it while expunging the “small”coefficients.

[0250] Obviously, such a mistake might be crucial, since the de-noisedsignal will have spurious, burst activity in finer resolution levels.Thus, the situation at hand is very contradictory. In order to preventpossible oversmoothing of the signal to be de-noised we have to set thethreshold ignoring noise samples having extremely large variance, butfrom the other hand, it might give rise to the serious estimation error.

[0251]FIG. 6 shows wavelet coefficients of the white noise with varianceσ²=1 (left-hand plot) and wavelet coefficients of the Fisher-Tippetnoise with variance σ²=0.25 (right-hand plot). Note that for both noisesthe dotted lines denotes the threshold computed according to (6.1) andt_(n={square root}) {square root over (2·log(n))}·σ. We see that all thewhite noise wavelet coefficients shown in the left-hand plot of thefigure are “captured” by the threshold guaranteeing us a noise-freeresult. In the same time, the most (almost all) Fisher-Tippet noisecoefficients are in the range to be shrunk, however there are fewcoefficients exceeding beyond the threshold. No doubts, that thesecoefficients can “spoil” all absolutely.

[0252] Whether the situation is hopeless? Not completely. In the nextchapter we will suppose a very simple and efficient solution for it.Before this, we note that the same problem exists in the imageprocessing applications considering the speckle noise reduction. It wasshown (H. H. Arsenault, G. April, 1976 and D. L. Goodman, 1976), thatthe coherent speckle noise can be modeled as a multiplicative noise,i.e. an image contaminated by the speckle noise is the original imagemultiplied by the white Gaussian noise. Consequently, after thelogarithmic transformation we have a sum of the logarithm of the originto be recovered and the Fisher-Tippet (“almost Gaussian”) noise. Thusthe goal for speckle reduction is equivalent to finding the bestestimate of the former by means of the wavelet shrinkage (i.e.,non-linear thresholding of the wavelet coefficients), for instance, asit was reported in [11-12].

[0253] 8. Non-linear Wavelet Analysis

[0254] The theory proposed by David Donoho and Iain Johnstone is basedon the assumption that the noise contaminating a function to berecovered is close to the Gaussian, white noise. As a result, theirthree-step de-noising algorithm is very sensitive to outliers. In ourspecific case it falls for the noise samples coming from the “long tail”of the Fisher-Tippet distribution and having relatively large variance.The problem is that the “outliers” are treated as local features, and,consequently, preserved.

[0255] To take around this problem Outlier resistant Wavelet transformswas proposed in [13]. The pyramid algorithm was adapted to the form inwhich outliers are captured into robust residuals at differentresolution levels. This algorithm was named Smoother-Cleaner WaveletTransform (SCWT). The idea is as follows. The DWT is a lineartransformation, and consequently its coefficients are very sensitive tooutliers. The classical Wavelet transform yields a sequence ofapproximations f_(k)(t) of the original function f(t). Theapproximations are nothing than the projections of the function into thesub-spaces spanned by the collection of scaling functions$\begin{matrix}{\varphi_{k,m} = {2^{\frac{k}{2}} \cdot {\varphi \left( {{2^{- k}t} - m} \right)}}} & (7.1)\end{matrix}$

[0256] In each sub-space each projection minimize the L₂ norm

∥f(t)−f _(k)(t)∥₂ ²  (7.2)

[0257] The problem is that the L₂ norm is extremely sensitive tooutliers. At this point we have to note that it is very crucial to theprojection method for estimation of the ultrasound pulse log-spectrum(see appropriate chapter). Let's remind that its estimate was obtainedas a projection of the RF-line log-spectrum on an appropriate subspace,i.e. the sub-space spanned by a collection the scaling functions (asit's given in (7.1)) with a prescribed index k. For the above- mentionedreason the outliers could significantly distort such an estimate. So thediscussion below is also relevant to the projection algorithm.

[0258] One way to reduce the sensitivity of the Wavelet transform tooutliers is to construct projections minimizing the L₁ norm (that iswell known to be quite resistant to them). It was shown that in thiscase outliers could be captured completely in few projections ontocomplement subspaces, i.e. sub-spaces spanned by the collection ofmother wavelets $\begin{matrix}{\psi_{k,m} = {2^{\frac{k}{2}} \cdot {\psi \left( {{2^{- k}t} - m} \right)}}} & (7.3)\end{matrix}$

[0259] In such a case it is easier to remove outliers from the L₁decomposition. Unfortunately, this evidently attractive decomposition isinefficient from the computational point of view. For this reason SCWTwas proposed. This Wavelet decomposition is very fast (computationalcomplexity O(n)), and it performs like the classical Wavelet transformfor Gaussian signals preventing outliers from leaking into the waveletcoefficients at coarser levels. The idea is to isolate waveletcoefficients probably ascribed to outliers in robust residuals at eachresolution level.

[0260] Let's S(k) denote wavelet coefficients at level k, and Ŝ(k) is arobust smooth of S(k) obtained by running median of length 5 (note thatthe length of the robust filter should be as small as possible toprovide minimal distortion of the underlying signal). The robustresiduals R(k) of the level k are computed as follows $\begin{matrix}{{R(k)} = \left\{ \begin{matrix}{0,} & {if} & {{{{S(k)} - {\overset{\Cap}{S}(k)}}} \leq \lambda_{k}} \\{{{{sign}\left( {{S(k)} - {\overset{\Cap}{S}(k)}} \right)} \cdot \left( {{{{S(k)} - {\overset{\Cap}{S}(k)}}} - \lambda_{k}} \right)},} & {if} & {{{{S(k)} - {\overset{\Cap}{S}(k)}}} \geq \lambda_{k}}\end{matrix} \right.} & (7.4)\end{matrix}$

[0261] The threshold λ_(k) is chosen so that most of the robustresiduals are zeroed. Subsequently the classical Wavelet decompositionstep resulting in detail and approximation wavelet coefficients S(k+1)are applied to [S(k)−R(k)]. Clearly that the above steps are repeatedwith S(k+1) and so on.

[0262] It was shown that in the Gaussian noise case, the proposedalgorithm produce essentially the same decomposition as the classicalDWT. When the signal to be analyzed is contaminated by a noise withoutliers and its patches, it is completely captures in waveletcoefficients of finer levels retaining the signal approximationundistorted.

[0263] 9. Estimation by Robust De-noising

[0264] Returning to the pulse log-spectrum estimation problem, let'sremind that it's equivalent to the problem of de-noising of the desiredsignal, when the noise to be rejected has the Fisher-Tippetdistribution. In previous sections it was shown that a sequence of thenoise could be viewed as a Gaussian white noise contaminated by“outliers” having extremely large variance. Generally, the residuals ofthe first decomposition level correspond to isolated outliers or pairsof outliers. In this case the simplest was to perform de-noising is todiscard the robust residuals and to use wavelet shrinkage standardalgorithm on the wavelet coefficients. The procedure is shown in thefollowing block-diagram.

[0265] Block-diagram No1: Robust Ultrasound Pulse Log-spectrumEstimation Algorithm Using Outlier Resistive De-noising

[0266] In other words, in order to estimate the pulse log-spectrum wecompute robust residual of a RF-line log-spectrum according to (7.4),then subtract it from the latter, and subsequently apply the classicalde-noising algorithm of David L. Donoho and Iain Johnstone to theobtained difference.

[0267] Note that if the robust median filter is not long (from 5 to 7points), then it has no any influence on the signal to be recovered.Consequently the above described subtraction and discard of residualsmust change the probability distribution of the noise. This situation isdepicted in FIG. 7. In its upper plot one can see a histogram of asimulated i.i.d. Fisher-Tippet sequence viewed at once with itstheoretical Gaussian approximation. In the lower plot of the picture onecan see the noise histogram after removing residuals, which werecomputed according to (7.4). Note the disappearance of the “long tail”indicating that all the “outliers” were removed. Additionally one cannote that after the above “outlier shrinkage” procedure the noisedistribution apparently resembles the Gaussian probability densityfunction with variance that is very close to be theoretically predictedvariance of the Gaussian approximation of the Fisher-Tippetdistribution, namely 0.25. It guarantees that the “modified” noise doesnot include samples, which could be treaded as features of the signal tobe recovered, and it is reasonable to expect that the three-stepde-noising algorithm will supply us desired results.

[0268] Before the demonstration of results we would like to add somenotes regarding choice of the threshold in (7.4) and wavelet. Inprincipal, we can compute the threshold value using theoretical value ofthe variance of the Gaussian approximation. Yet, from the practicalpoint of view it was observed that one could obtain quite acceptableresults choosing the threshold in such a way that 90-95% of thedifferences [S(k)−Ŝ(k)] are under it, and as a result are annihilated.Regarding the DWT filter choice it was observer that nearly symmetricfilters of I. Daubechies perform very well.

[0269] Now let's show results of our computer simulations. In FIG. 8 thenecessity of using the smoothing-shrinkage step turning thereconstruction to be insensitive to “large-variance” samples of theFisher-Tippet noise is demonstrated. The left-hand plot shoes theestimated pulse log-spectrum, when the estimation was performed usingthe classical three-step de-noising algorithm applied to the RF-linelog-spectrum. One can see that indeed the noise “outliers” are treatedas local features of the signal to be recovered, and consequentlypreserved in the reconstruction in the form of spurious details in finerresolution levels. In the right-hand plot one can see the samereconstruction preceded by the outliers smoothing-shrinkage step.Obviously there are no evident signs of influence the outliers on theestimate quality. From the quantitative point of view, the normalizedmean-squared error in the first case is 4.736%, whereas in the secondcase 3.982%. Thus, one can conclude that the de-noising technique is apowerful tool for the pulse spectrum estimation, when it is preceded bythe appropriate treatment of the noise “large-variance” samples.

[0270] The proposed algorithm uses the decomposition of RF-linelog-spectra in the dyadic Wavelet basis. However, it should be notedthat the performance of the pulse estimation algorithm might besignificantly improved by using another kinds of transformations, likeWavelet Packet, Trigonometric bases, Shift-invariant libraries and etc.,i.e., as it was described in [1-4]. So, many interesting questionsremained out of scope of the paper. What results would be obtained byapplying the technique in the “best” basis domain? WhetherTranslation-Invariant De-Noising can improve the pulse reconstructionquality? . . . The work is not completed and the problem deservesfurther study.

[0271] 10. References

[0272] [1] “De-noising by Soft-Thresholding”, David L. Donoho, 1992

[0273] [2] “Adapting to Unknown Smoothness via Wavelet Shrinkage” DavidL. Donoho and Iain Johnstone, 1993

[0274] [3] “Ideal De-noising in an Orthonormal Basis Chosen from aLibrary of Bases”, David L. Donoho and Iain Johnstone, 1994

[0275] [4] “Translation-Invariant De-noising”, David L. Donoho and R. R.Coifman, 1995

[0276] [5] “Wavelet shrinkage and Wavelet Vaguelette Decomposition: ATen-Minute Tour”, David L. Donoho, 1993

[0277] [6] “Nonlinear Wavelet Methods for Recovering Signals, Images,and Densities from indirect and noisy data”, David L. Donoho, 1993

[0278] [7] “Wavelets on the intervals and fast wavelet transforms”, ACohen, I. Daubechies, P. Vial, 1993

[0279] [8] “Properties of speckle integrated with a finite aperture andlogarithmically transformed”, H. H. Arsenault, G. April, Journal ofOptical Society of America, 66, November, 1976

[0280] [9] “Some fundamental properties of speckle”, D. L. Goodman,Journal of Optical Society of America, 66, November, 1976

[0281] [10] “Nonlinear processing of a shift invariant DWT for noisereduction”, M. Lang, H. Guo, J. E. Odegard, and C. S. Burrus, SPIE, Vol.2491, 1995

[0282] [11] “Simultaneous speckle reduction and data compression usingbest wavelet packet bases with application to SAR based ATD/R”, D. Wei,J. E. Odegard, M. Lang, and C. S. Burrus, SPIE, Vol. 2491, 1995

[0283] [12] “Speckle reduction via Wavelet Shrinkage with application toSAR based ATD/R”, H. Guo, J. E. Odegard, M. Lang, R. A Gopinath, I. W.Selesnick, and C. S. Burrus, SPIE, Vol. 2303, 1994

[0284] [13] “De-noising and robust non-linear Wavelet analysis”, AndrewG. Bruce, David L. Donoho, Hong-Ye Gao, R. Douglas Martin, SPIE, Vol.2242, 1994

[0285] [14] “Nonparametric estimation of ultrasonic pulses using the FWTand determination of the reflectivity distribution using the approximateinverse technique”, D. Adam, O. Michailovich, 2000, to be appeared.

What is claimed is:
 1. A method of imaging a target, comprising thesteps of: (a) acquiring an echo sequence image of the target; (b)computing a log spectrum of at least a portion of said echo sequenceimage; (c) computing a low-resolution wavelet projection of said logspectrum; (d) estimating a point spread function from saidlow-resolution wavelet projection; and (e) deconvolving said at leastportion of said echo sequence image with said point spread function. 2.The method of claim 1, further comprising the step of: (f) partitioningthe echo sequence image into a plurality of segments; said computing ofsaid log spectrum, said computing of said wavelet projection, saidestimating and said deconvolving then being effected separately for eachsaid segment.
 3. The method of claim 2, wherein said segments aredisjoint.
 4. The method of claim 2, wherein at least two of saidsegments at least partly overlap.
 5. The method of claim 1, wherein saidwavelet projection is based on a wavelet selected from the groupconsisting of Coiflet wavelets, minimum-phase Daubechies wavelets,symmetric Daubechies wavelets and biorthogonal wavelets.
 6. The methodof claim 1, wherein said computing of said wavelet projection includeseffecting an outlier-resistant wavelet transform of said log spectrum.7. The method of claim 1, wherein said estimating of said point spreadfunction includes estimating a frequency domain phase of said pointspread function.
 8. The method of claim 7, wherein said phase estimatingassumes that said point spread function is a minimum phase sequence. 9.The method of claim 1, wherein said deconvolving is effected using anapproximate inverse.
 10. A method of imaging a target, comprising thesteps of: (a) acquiring an echo sequence image of the target; (b)computing a log spectrum of at least a portion of said echo sequenceimage; (c) effecting an outlier-resistant wavelet transform of said logspectrum, thereby producing a plurality of wavelet coefficients; (d)soft-thresholding said wavelet coefficients; (e) applying an inversewavelet transform to said soft-thresholded wavelet coefficients toobtain a point spread function log spectrum; (f) estimating a pointspread function from said point spread function log spectrum; and (g)deconvolving said at least portion of said echo sequence image with saidpoint spread function.
 11. The method of claim 10, further comprisingthe step of: (h) partitioning the echo sequence image into a pluralityof segments; said computing of said log spectrum, said effecting of saidwavelet transform, said soft-thresholding, said applying, saidestimating and said deconvolving then being effected separately for eachsaid segment.
 12. The method of claim 11, wherein said segments aredisjoint.
 13. The method of claim 11, wherein at least two of saidsegments at least partly overlap.
 14. The method of claim 10, whereinsaid outlier-resistant wavelet transform is based on minimizing an L₁norm.
 15. The method of claim 10, wherein said outlier-resistant wavelettransform is a Smoother-Cleaner Wavelet Transform.
 16. The method ofclaim 10, wherein said outlier-resistant wavelet transform includes aplurality of resolution levels, with robust residuals being removed inonly a portion of said resolution levels.
 17. The method of claim 16,wherein said robust residuals are removed only in a first saidresolution level.
 18. The method of claim 10, wherein said estimating ofsaid point spread function includes estimating a frequency domain phaseof said point spread function.
 19. The method of claim 18, wherein saidphase estimating assumes that said point spread function is a minimumphase sequence.
 20. The method of claim 10, wherein said deconvolving iseffected using an approximate inverse.
 21. An apparatus for imaging atarget, comprising: (a) a transducer for acquiring an echo sequenceimage of the target; and (b) a processor for: (i) computing a logspectrum of at least a portion of said echo sequence image, (ii)computing a low-resolution wavelet projection of said log spectrum,(iii) estimating a point spread function from said low-resolutionwavelet projection, and (iv) deconvolving said at least portion of saidecho sequence image with said point spread function.
 22. An apparatusfor imaging a target, comprising: (a) a transducer for acquiring an echosequence image of the target; and (b) a processor for: (i) computing alog spectrum of at least a portion of said echo sequence image, (ii)effecting an outlier-resistant wavelet transform of said at leastportion of said log spectrum, thereby producing a plurality of waveletcoefficients, (iii) soft-thresholding said wavelet coefficients, (iv)applying an inverse wavelet transform to said soft-thresholded waveletcoefficients to obtain a point spread function log spectrum, (v)estimating a point spread function from said point spread function logspectrum, and (vi) deconvolving said at least portion of said echosequence image with said point spread function.